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ABSTRACT 


Phase -shift  keyed  (PSK)  signal  modulation  methods  are 
coming  into  increasing  use  in  modern  communications.  This 
thesis  describes  the  implementation  of  three  methods  of 
computing  the  cyclic  spectrum  to  determine  the  presence  of  PSK 
signals.  The  Strip  Spectral  Correlation  Algorithm  (SSCA)  and 
the  Fast  Fourier  Transform  (FFT)  Accumulation  Method  (FAM) 
both  estimate  the  cyclic  spectral  plane.  The  Sub- FFT 
Accumulation  Method  (SUBFAM)  program  computes  the  Spectral 
Correlation  Function  (SCF)  for  a  set  of  possible  spectral 
frequencies  for  a  single  cycle  frequency.  The  results  of 
algorithm  performance  are  presented  along  with  a  discussion  of 
promising  areas  for  perfor.aance  enhancement  and  automation  of 
signal  detection  and  classification. 
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I .  INTRODUCTION 


A.  BACKGROUND 

The  field  of  cyclic  spectral  analysis  is  growing  in 
importance  as  non- stationary  modulation  schemes  proliferate  in 
modern  communications.  Traditionally,  spectral  analysis 
methods  have  assumed  that  signals  of  interest  were  stationary 
or  had  non- time -varying  statistical  parameters.  This 
asstimption  is  untrue  for  the  cyclostationary  signals  created 
when  using  pulse  or  carrier  amplitude  modulation,  phase  or 
frequency  carrier  modulation  or  digital  pulse  or  carrier 
modulation.  Gardner  [Ref.  l:pp  419-453]  discusses  many  of 
these  signal  types.  Spread  spectrum  modulation  techniques 
such  as  direct  sequence  PSK  and  frequency -hopped  FSK  also 
exhibit  temporally  varying  spectral  features  [Ref.  l:pp  453- 
457]  .  Reference  2  provides  a  tutorial  on  cyclic  spectral 
analysis . 

Cyclic  spectral  analysis  methods  take  advantage  of 
spectral  variation  over  time  by  correlating  spectral 
estimations  at  discrete  time  intervals  to  locate  signals  that 
may  not  otherwise  be  identified.  A  direct  sequence  signal 
with  a  very  wide  bandwidth  may  be  indistinguishable  from  a 
noisy  background  in  the  frequency  spectral  estimate  of 
Figure  1  [Ref.  3:pp  2-9  -  2-10].  But  this  signal  becomes 
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clearly  identifiable  in  the  cyclic  spectral  estimate  of 
Figure  2 . 

One  cyclic  spectral  technique  is  the  frequency  smoothed 
cyclic  periodogram  method  (FSM)  [Ref.  l:pg  464].  This  has 
been  implemented  in  the  Cyclic  Spectral  Analysis  Software 
Package  [Ref.  4]  as  the  program  sxaf.  Hereafter,  this 
software  will  be  referred  to  as  the  SSPI  package  or  the  FSM. 
This  method  was  used  to  generate  Figures  1  and  2 . 
Unfortunately,  the  FSM  is  less  efficient  than  time  smoothing 
methods  for  general  spectral  estimation  [Ref.  5:pg  38]. 

B.  THESIS  GOALS 

The  purpose  of  this  thesis  is  to  implement  two  time 
smoothing  algorithms,  the  FFT  Accumulation  Method  (FAM) 
[Ref.  5:pp  44-47]  and  the  Strip  Spectral  Correlation  Algorithm 
(SSCA)  [Ref.  5:pp  47-48].  The  FAM  is  less  computationally 
complex  than  the  FSM  and  the  SSCA  less  so  than  the  FAM.  Both 
methods  have  been  written  to  be  compatible  with  the  SSPI 
package.  By  integrating  FAM  and  SSCA  into  the  SSPI  package 
future  research  may  more  easily  be  done  in  algorithm 
performance  comparison  and  enhancement . 

Both  the  FAM  and  SSCA  accept  input  data  generated  by  the 
SSPI  utility  program  pam.  They  also  generate  results 
compatible  with  the  SSPI  plotting  utility  plot_sxaf.  Each 
program  may  optionally  generate  output  in  a  format  compatible 
with  the  MATLiAB  [Ref.  6]  package  which  is  in  common  use  at  the 


Naval  Postgraduate  School.  Input  and  output  data  are  in  ASCII 
file  formats. 

Throughout  this  thesis  the  same  BPSK  signal  data  is  used 
to  illustrate  the  implementations  of  FAM  and  SSCA.  A 
representative  signal  sample  is  plotted  in  Figure  3  .  No  noise 
was  added  to  this  signal.  By  estimating  the  frequency 
spectrum  at  different  points  in  time,  it  is  evident  in  Figure 
4  that  the  spectral  features  do  indeed  vary  temporally. 

Chapter  II  describes  the  implementation  of  the  FAM. 
Chapter  III  describes  the  implementation  of  the  SSCA.  Chapter 
IV  describes  a  utility  to  estimate  cycle  frequencies  for  a 
given  baseband  frequency.  The  Sub-FFT  Accumulation  Method 
(SUBFAM)  program  [Ref.  7]  gives  a  quick  result  for  a  specific 
subset  of  the  spectral  plane.  Chapter  V  discusses  specific 
algorithm  performance  comparisons.  Finally,  Chapter  VI 
presents  conclusions  and  recommendations  for  future  areas  of 
research. 
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Figiirs  3  Typical  BPSK  Signal  Data 


Figure  4  Tine  Sequence  Speccrum  of  Typical  BPSK  Signal 
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II.  FFT  ACCUMULATION  METHOD 


A.  THEORY 

The  FFT  Accumulation  Method  (FAM)  [Ref.  5:pp  44-47]  is  a 
Fourier  transform  of  correlation  products  between  spectral 
components  smoothed  over  time.  Periodicities  in  the  spectral 
components  then  become  detectable.  The  cyclic  spectral  plane 
as  shown  in  Figure  5  ranges  in  frequency  from  -tj2  to  +fs/2, 
and  in  cycle  frequency  from  -f,  to  +f5  where  f^  is  the  sampling 
frequency.  For  each  unique  cell  in  Figure  5,  the  FAM 

cross- spectral  estimate  is  defined  to  be  [Ref.  5:pg  44]: 

P-l  -i2Kza 

r;(rL,f,)g^(n-r)  e  ^  (l) 

r«o 


I 

The  FAM  auto- spectral  estimate  is  defined  to  be 
[Ref .  5 :pg  44] : 


^XX-r 


p-l 

{nL,  fj)  =  -ffc)  Xj(rL,  f,)  g^[n-r)  e 


r^Q 


-i2nrq 
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where:  j  is  the  average  frequency  bin  (k;+t’)/2 

q  is  the  difference  or  bandwidth  frequency 
g,  (n-r)  represents  the  Hamming  window. 

The  FAM  was  implemented  by  forming  an  array  from  x(kT) 
(0  <  k  s  N-1)  with  rows  which  are  N'  points  long  from  the 
input  sample  data.  The  starting  point  of  each  succeeding  row 
is  offset  from  the  previous  rows  starting  position  by  L 
samples.  A  Hamming  window  is  applied  across  each  row  which  is 
then  Fast  Fourier  transformed  and  downconverted  to  baseband. 
The  result  at  this  point  is  a  two-dimensional  array  with 
columns  representing  constant  frequencies.  Each  column  was 
point-wise  multiplied  in  turn  with  the  conjugate  of  the 
columns  resulting  from  processing  y(kT)  (0  s  k  s  N-1) .  Each 
resultant  product  vector  was  Fast  Fourier  transformed  and  the 
low  frequency  half  placed  into  the  final  cyclic  spectral  plane 
at  the  appropriate  location  in  the  cyclic  spectral  plane. 
Figure  6  shows  a  block  diagram  for  the  FAM  cross -spectral 
estimate . 

The  following  subsections  in  this  chapter  discuss  in 
detail  what  has  been  described  in  the  previous  paragraph.  The 
cycle  frequency  resolution  of  FAM  is  Aa  =  fJPL  [Ref.  5:pg  44] 
where  f^  is  the  original  data  sampling  rate.  The  frequency 
resolution  of  FAM  is  Af  =  fjN’  [Ref.  5:pg  44]  .  The  time- 
frequency  resolution  product  is  AtAf  =  PL/N'  [Ref.  5:pg  45]. 
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The  command  line  format  for  calling  the  FAM  program  is 
provided  in  Appendix  A.  The  FAM  source  code  listing  is 
provided  in  Appendix  B. 

B .  IMPLEMENTATION 

1.  Input  Channelization 

The  input  sample  data  is  formed  into  a  two-dimensional 
array.  The  array  row  length  is  equal  to  the  number  of  input 
channels  N'  .  For  a  given  number  of  input  sample  points  N,  a 
row  size  of  N’  ,  and  a  chosen  offset  L,  there  are  P  =  (N- 
N' ) /L  =  N/L  rows  formed.  The  choice  of  N'  must  take  into 
consideration  that  ideally  the  t  ime  -  f r equency  resolution 
product  must  be  much  greater  than  one  [Ref.  5:pg  40].  N’ 
should  also  be  a  power  of  2  to  avoid  truncation  or  zero¬ 
padding  in  the  FFT  routines.  L  should  be  chosen  to  be  less 
than  or  equal  to  N'  /4 . 

The  completely  filled  array  is  P  rows  by  N'  columns. 
Figure  7  shows  how  a  small  array  is  filled  from  a  discretely 
sampled  signal  x(kT)  when  N'=16,  P=8  and  L=4.  The  number 
inside  each  cell  represents  the  value  of  k  used  to  index  on 
x(kT)  to  fill  that  location  in  the  array. 

Figure  8  shows  the  magnitude  of  the  original  data  for 
the  example  BPSK  signal  organized  into  the  P  by  N'  array  where 
N=512,  N'=32,  L=4  and  P=128.  The  input  data  is  assumed  to  be 
complex  with  a  real  and  imaginary  component  to  each  sample 
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point.  Figure  9  shows  a  single  row  of  the  array.  The  phase 
changes  of  the  BPSK  are  evident . 

2 .  Windowing 

A  Haituning  window  [Ref  7:pg  467]  is  applied  to  each  row 
of  the  array.  The  equation  for  the  Haniming  window  is: 

w(r)  =0.54-0. 46cos(-^il^)  ,  0  ^  r  ^  N^-1  (3) 

N^-1 

A  32 -point  Hamming  window  is  plotted  in  Figure  10.  It 
is  applied  to  both  the  real  and  imaginary  parts  of  the  complex 
example  array.  The  magnitude  of  the  resultant  array  is  shown 
in  Figure  11. 

3 .  First  PPT 

Each  row  of  the  windowed  data  array  is  Fast  Fourier 
transformed  to  reveal  the  first  spectral  components.  The 
resultant  array  is  still  indexed  P  rows  by  N'  columns  but  now 
the  column  index  relates  to  a  specific  bin  of  spectral 
frequencies.  Figure  12  illustrates  this  relationship. 
Figure  13  shows  the  results  of  FFTing  the  BPSK  example. 

4 .  Downconversion 

Each  row  of  spectral  components  is  downconverted  to 
baseband  through  multiplication  with  the  complex  exponential, 


-iZnkmL 


where:  m  is  the  row  index,  0  sms  P-i 

k  is  the  column  index,  0  s  k  s  iV' - 1 

The  magnitude  of  the  exponential  is  unity  over  the 
array  but  the  phase  shows  considerable  variation.  Figure  14 
shows  the  phase  of  the  exponential  over  the  P  by  N'  array. 
Figure  15  shows  the  phase  for  one  representative  row.  The 
magnitude  of  the  array  remains  unchanged  from  Figure  13 . 

5.  Multiplication 

For  each  cell  in  the  cyclic  spectral  plane,  one  column 
of  the  array  is  multiplied  with  the  conjugate  of  another 
column.  The  separation  of  the  columns  is  determined  by  the 
desired  frequency  separation  or  cycle  frequency  {cy=f|j-f,)  .  The 
mid-point  between  the  columns  is  the  frequency  bin  of  interest 
(f=  (f^+f,) /2)  .  Figure  16  shows  two  columns  to  be  conjugate 
multiplied  for  use  in  filling  a  specific  f,a  cell.  Figure  17 
illustrates  the  f,a  bin  values  in  the  cyclic  spectral  plane 
for  W' =4 .  Figure  18  shows  the  corresponding  two  column 
indices  used  to  form  the  product  vector  which  is  used  to 
produce  the  cells  of  Figure  17. 

6 .  Second  FFT 

The  product  from  the  previous  multiplication  is  FFT'd 
to  yield  a  P-point  result.  Only  the  middle  of  the  resulting 
spectrum  is  retained  and  stored  into  the  designated  f.a  cell. 
The  upper  and  lower  ends  are  undesireable  because  of  increased 
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estimate  variation  at  the  channel  pair  ends  [Ref.  5:pg  45]. 
Figure  19  illustrates  which  part  of  the  estimate  is  retained 
and  placed  into  the  cell.  Figure  20  shows  the  final  result  of 
the  FAM  auto- spectral  process  on  the  example  BPSK  signal. 

7 .  Data  Reduction 

The  FAM  program  typically  generates  large  output  data 
files.  For  convenience,  an  option  may  be  chosen  to  reduce  the 
amount  of  output.  By  comparison  sorting  for  the  largest  a 
value  in  an  f,0!  cell,  the  number  of  a  slices  output  is  reduced 
from  N'V/2  to  N' .  Overall  FAM  program  execution  time 
increases  accordingly  to  accomplish  the  sorts. 

Figure  21  illustrates  the  output  full  cyclic  spectral 
plane  storage  array  before  sorting.  Figure  22  shows  the 
output  array  after  data  reduction  has  been  completed.  Figures 
23  and  24  plot  the  resulting  spectral  half -planes  without  and 
with  data  reduction  respectively. 

Since  all  cells  have  an  equal  number  of  a  values,  all 
sorts  are  of  equal  complexity.  Further  work  on  data  reduction 
would  require  selection  of  a  largest  a  value  from  widely 
varying  numbers  of  a  values  depending  on  the  choice  of  N,  N' 
and  L. 

C .  PERFORMANCE 

An  established  method  of  evaluating  the  complexity  of  an 
algorithm  is  to  determine  the  number  of  floating  point 


operations  that  must  be  performed.  For  this  estimate  it  is 
assumed  that  an  auto-spectral  estimate  is  being  computed  and 
the  data  is  complex  in  every  step.  It  is  also  assumed  that 
the  Hamming  window  coefficients  and  the  downconversion 
coefficients  have  been  previously  computed  and  stored  for 
later  use.  Each  N-point  FFT  requires  {N/2)*log2N  complex 
floating  point  multiplies  [Ref.  8:pg  506]  or  4* (N/2 ) *log2N 
real  floating  point  multiplies.  The  cost  of  any  output  data 
reduction  is  not  considered  here. 


Applying  the  window 

First  FFT . 

Downconversion . 

Multiplication . 

Second  FFT . 


2*P*W' 

2*p*N’  *log,W' 
4*P*W' 

4*W'  *N'  *P 
2*N’  *N’  *P*log2P 


Total:  (6+4*W' )  *P*W'  +  (2*P*W' )  *  dog^W'  +  W'^log^P)  (4) 

Since  P-N/L  (5) 

Total:  (6+4*W' )  *NW'/L  +  (2*Ni^' /L)  ( logjW'  +  N'logjN/L)  (6) 
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Figure  7  Layout  of  A  Sample  Input  Array  Showing  Input 
Sample  Storage  for  N' =lo ,  N=48,  L=4  and  ?=8. 
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igure  11  Example  3PSK  Signal  Data  After  Window, 
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Figure  19  Retained  Section  of  the  Second  FFT  Results 
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Figure  20 
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III.  STRIP  SPECTRAL  CORRELATION  ALGORITHM 


A.  THEORY 

The  Strip  Spectral  Correlation  Algorithm  (SSCA) 
[Ref.  5:pp  47-48]  is  a  Fourier  transform  of  correlation 
products  between  spectral  and  temporal  components  smoothed 
over  time.  Periodicities  in  the  spectral  components  then 
become  detectable.  The  cyclic  spectral  plane  as  shown  in 
Figure  25  ranges  in  frequency  from  -tj2  to  +fj2,  and  in  cycle 
frequency  from  -f,  to  +f,  where  f^  is  the  sampling  frequency. 
For  each  unique  (fj,a<,)  strip  in  Figure  25,  the  SSCA  cross - 
spectral  estimate  is  defined  to  be  [Ref.  5:pg  47]: 


-iZitrq 

n 


(7) 


r«o 


The  SSCA  auto-spectral  estimate  is  defined  to  be 
[Ref .  5 :pg  47]  : 


^XXj. 


<"-T- 


gAa 


73-1 


Xj.{r,  f^)  X'  (r)  g{n~r)  e 


-i2nzq 

n 


(8) 


r-O 


where:  k  is  a  multiple  of  the  frequency  resolution, 

-N'  /2  s  k  s  (N' /2)  -1 

q  is  a  multiple  of  the  cycle  frequency  resolution, 
-N'  s  q  s  (17'  -1) 

g{n-r)  represents  the  Hamming  window. 

The  SSCA  was  implemented  by  forming  an  array  from  x(kT) 
(0  s  k  £  N-1)  with  rows  which  are  N'  points  long  from  the 
input  sample  data.  The  starting  point  of  each  succeeding  row 
is  offset  from  the  previous  rows  starting  position  by  L 
samples.  A  Hamming  window  is  applied  across  each  row  which  is 
then  Fast  Fourier  transformed  and  downconverted  to  baseband. 
The  result  at  this  point  is  a  two-dimensional  array  with 
columns  representing  constant  frequencies.  Each  row  is 

transposed  and  replicated  L  times  for  a  total  of  PL  columns. 
Each  column  is  then  multiplied  by  a  sample  y(kT)  (0  s  k  £  (PL- 
1)  )  .  Each  resultant  row  of  the  array  is  Fast  Fourier 

transformed  and  placed  into  a  strip  of  the  final  cyclic 
spectral  plane  at  the  appropriate  location.  Figure  26  shows 
a  block  diagram  for  the  SSCA  cross -spectral  estimate. 

The  following  subsections  in  this  chapter  discuss  in 

detail  what  has  been  described  in  the  previous  paragraph.  The 
cycle  frequency  resolution  of  SSCA  is  Aa  =  fj/N  [Ref.  5:pg  48] 
where  f,  is  the  original  data  sampling  rate  and  N  is  the 

number  of  input  samples  used.  The  frequency  resolution  of 
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SSCA  is  equal  to  the  sampling  rate  divided  by  the  number  of 
channels  N'  ,  Af  =  tjN'  [Ref.  5:pp  42,48]  .  The  time- frequency 
resolution  product  is  AtAf  =  N/i\^'  [Ref.  5:pg  48]  . 

The  command  line  format  for  calling  the  SSCA  program  is 
provided  in  Appendix  C.  The  SSCA  source  code  listing  is 
provided  in  Appendix  D. 

B .  IMPLEMENTATION 

1.  Input  Channelization 

The  input  sample  data  is  formed  into  a  two-dimensional 
array .  The  array  row  length  is  equal  to  the  number  of  input 
channels  N‘  .  For  a  given  number  of  input  sample  points  N,  a 
row  size  of  N'  ,  and  a  chosen  offset  L,  there  are  P  =  (N- 
N' )  JL  -•  N/L  rows  formed.  The  choice  of  N’  must  take  into 
consideration  that  ideally  the  time- frequency  resolution 
product  must  be  much  greater  than  one  [Ref.  5:pg  40].  N’ 
should  also  be  a  power  of  2  to  avoid  truncation  or  zero- 
padding  in  the  FFT  routines.  L  should  be  chosen  to  be  less 
than  or  equal  to  N' /4 . 

The  completely  filled  array  is  P  rows  by  N’  columns. 
Figure  27  shows  how  a  small  array  is  filled  from  a  discretely 
sampled  signal  x(kT)  when  N'=16,  P=8  and  L=4.  The  number 
inside  each  cell  represents  the  value  of  k  used  to  index  on 
x(kT)  to  fill  that  location  in  the  array. 

Figure  28  shows  the  magnitude  of  the  original  data  for 
the  example  BPSK  signal  organized  into  the  P  by  N'  array  where 
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N=512,  N'=32,  L=4  and  P=128.  The  input  data  is  assumed  to  be 
complex  with  a  real  and  imaginary  component  to  each  sample 
point.  Figure  29  shows  a  single  row  of  the  array.  The  phase 
changes  of  the  BPSK  are  evident. 

2 .  Windowing 

A  Hamming  window  [Ref  7:pg  467]  is  applied  to  each  row 
of  the  array.  The  equation  for  the  Hamming  window  is: 

w{  r)  =0.54-0. 46  cos(  ),  0  ^  r  ^  N'-1  (9) 

N‘-1 

A  32 -point  Hamming  window  is  plotted  in  Figure  30.  It 
is  applied  to  both  the  real  and  imaginary  parts  of  the  complex 
example  array.  The  magnitude  of  the  resultant  array  is  shown 
in  Figure  31. 

3 .  First  FPT 

Each  row  of  the  windowed  data  array  is  Fast  Fourier 
transformed  to  reveal  the  first  spectral  components.  The 
resultant  array  is  still  indexed  P  rows  by  W'  columns  but  now 
the  column  index  relates  to  a  specific  bin  of  spectral 
frequencies.  Figure  32  illustrates  this  relationship. 
Figure  33  shows  the  results  of  FFTing  the  BPSK  example. 

4.  Oownconverslon 

Each  row  of  spectral  components  is  downconverted  to 
baseband  through  multiplication  with  the  complex  exponential, 
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where:  m  is  the  row  index,  0  sms  P-1 

k  is  the  column  index,  0  s  k  s  W' - 1 

The  magnitude  of  the  exponential  is  unity  over  the 
array  but  the  phase  shows  considerable  variation.  Figure  34 
shows  the  phase  of  the  exponential  over  the  P  by  N'  array. 
Figure  35  shows  the  phase  for  one  representative  row.  The 
magnitude  of  the  array  remains  unchanged  from  Figure  33. 

5.  Replication 

Each  row  is  copied  into  one  column  of  an  empty  N'  by 
PL  array.  It  is  then  replicated  in  the  L-l  adjacent  columns. 
Figure  36  illustrates  this  process. 

6.  Multiplication 

Each  column  of  the  array  is  pointwise  multiplied  with 
the  conjugate  of  a  sample  value  y(kT) .  There  are  PL«N  columns 
and  PL  samples  from  y(kT).  Figure  37  illustrates  the 
conjugate  multiplication  process. 

7 .  Second  FFT 

Each  row  from  the  previous  multiplication  is  Fast 
Fourier  transformed  to  yield  a  PL-point  result.  Each 
resulting  vector  is  stored  into  a  strip  of  the  cyclic  spectral 


plane.  Figure  38  shows  the  final  result  of  the  SSCA  auto- 
spectral  process  on  the  example  BPSK  signal. 

8 .  Data  Reduction 

The  SSCA  Program  typically  generates  large  output  data 
files.  For  convenience,  an  option  may  be  chosen  to  reduce  the 
amount  of  output.  By  comparison  sorting  for  the  largest  a 
value  in  an  f,Q;  cell,  the  number  of  a  slices  is  reduced  from 
N  to  N/L.  Overall  SSCA  execution  time  increases  accordingly 
to  accomplish  the  searches. 

Figure  39  illustrates  the  output  full  cyclic  spectral 
plane  storage  array  before  sorting.  Figure  40  shows  the 
output  array  after  data  reduction  has  been  completed.  Figures 
41  and  42  plot  the  resulting  spectral  half -planes  without  and 
with  data  reduction  respectively. 

Since  all  cells  have  an  equal  number  of  a  values,  all 
sorts  are  of  equal  complexity.  Further  work  on  data  reduction 
would  require  selection  of  a  largest  o;  value  from  widely 
varying  numbers  of  a  values  depending  on  the  choice  of  N,  N' 
and  L. 

C .  PERFORMANCE 

An  established  method  of  evaluating  the  complexity  of  an 
algorithm  is  to  determine  the  number  of  floating  point 
operations  that  must  be  performed.  For  this  estimate  i'  is 
assumed  that  an  auto-spectral  estimate  is  being  computed  and 
the  data  is  complex  in  every  step.  It  is  also  assumed  that 
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the  Hamming  window  coefficients  and  the  downconversion 
coefficients  have  been  previously  computed  and  stored  for 
later  use.  Each  N-point  FFT  requires  (N/2)*log;N  complex 
floating  point  multiplies  [Ref.  8:pg  506]  or  2*N*log2N  real 
floating  point  multiplies.  The  cost  of  any  output  data 
reduction  is  not  considered  here. 


Applying  the  window . 2*P*N' 

First  FFT . 2*P*N'  *log.N' 

Downconversion . 4  *  P*  W' 

Multiplication . A*N*N' 

Second  FFT . 2*N*W'*log2N 


Total : 

2*N'  * 

{e*P  +  4*N  +  (2*P 

+  2*N)*log2N) 

(10) 

since 

P-N/L 

(11) 

Total : 

2*N'  * 

(  (6*N/L  -t-  4*N)  +  1 

(2*N/L  +  2*N)*log2N) 

(12) 
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Figure  26  Block  Diagram  of  SSCA  Cross -Spectral  Estimate 
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Figxirs  27  Layout  of  A  Sample  Input  Array  Showing  Input 
Sample  Storage  for  W-lS,  N-4a,  L-4  and  P»a . 
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Figure  31  Example  BPSK  Signal  Data  After  Windowing 
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Figure  34  Phase  of  che 
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Figure  37  Multiplica 


Figure  39  3SCA  Storage  Array  Before  Data  Reduction 
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IV.  SUB-FFT  ACCLMDLATION  METHOD 


A.  THEORY 

In  Che  course  of  this  thesis  work  a  need  arose  for  a 
flexible  and  customized  Spectral  Correlation  Function  (SCF) 
[Ref.  7]  to  determine  the  presence  of  phase- shift  keyed 
signals  in  experimental  data.  The  SCF  [Ref.  7]  is 
represented: 


N-l 


S^(f)  =5^  x(r)  y(r)  e 


r-  2 


OSF^OSFye 


(13 


1  =  0 


where : 

N  is  the  number  of  data  samples  used 

f,  is  the  intermediate  modulation  frequency 

f,  is  the  sampling  frequency 

OSF  is  a  one-sided  bandpass  filter 

The  program  written  to  implement  this  algorithm  is  named 
the  Sub-FFT  Accumulation  Method  (SUBFAM)  program  to 
distinguish  it  from  the  more  general  purpose  FAM  program.  The 
command  line  format  for  calling  the  SUBFAM  program  is  provided 
in  Appendix  E.  The  SUBFAM  source  code  is  provided  in 
Appendix  F. 
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B .  IMPLEMENTATION 


Figure  43  illustrates  the  operations  performed  in  the 
SUBFAM  program.  A  set  of  N  points  is  obtained  from  within  a 
file  of  signal  data  values.  Each  set  of  N  points  is  processed 
through  a  one-sided  bandpass  filter  to  remove  either  all 
positive  or  negative  frequencies  as  desired.  Downconversion 
of  each  file  to  baseband  from  the  original  sampling  and 
transmission  frequency  bands  follows.  The  results  are  then 
correlated  through  a  complex  conjugate  multiplication.  A 
final  N-point  Fast  Fourier  Transform  yields  the  spectral 
correlation  result. 

C .  PERFOKMANCE 

The  SUBFAM  program  performed  as  required.  Figure  44 
illustrates  the  results  of  processing  test  data  containing 
phase- shift -keyed  signal  samples.  The  peaks  at  the  chip 
frequencies  are  apparent.  Figure  45  shows  the  results  of 
correlating  test  data  with  data  containing  only  white  noise. 
The  lack  of  correlation  between  spectral  features  yields 
results  which  are  four  orders  of  magnitude  less  than  in  the 
successful  signal  detection  of  Figure  44. 
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Figure  44  SUBFAM  Program  Result  With  Signal,  N-409S 
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V.  ALGORITHM  PERFORMANCE 

The  complexity  of  the  FAM  and  SSCA  algorithms  is  estimated 
in  Chapter  II  Section  C  and  Chapter  III  Section  C 
respectively.  Both  require  on  order  of  Nlog2N  floating  point 
multiplies  where  N  is  the  original  number  of  sample  points 
used.  Brown  and  Loomis  [Ref.  9]  explore  the  complexity  of  the 
FAM,  SSCA  and  FSM  algorithms  in  some  detail.  Their  results 
are  consistent  with  an  order  NlogjN  multiplies  for  FAM  and 
SSCA. 

The  cyclic  spectral  plane  is  symmetric  about  the  f=0  and 
the  Qi=0  axes  for  real  signals.  For  the  most  efficient 
performance  it  is  only  necessary  to  compute  one  quadrant  of 
the  cyclic  spectral  plane.  To  compare  the  performance  of  FAM 
and  SSCA  with  FSM,  complexity  figures  derived  from  Reference 
9  for  quadrant  complexity  will  be  used.  The  three  algorithms 
require  the  following  numbers  of  floating  point  multiplies: 

FAM  Complexity  [Ref.  9:pg  18] 

2Ni\r^log2(  — )  +  SNlog^CiyO  4NW''  -  20N  (14) 

n' 
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SS'CA  Complexity  [Ref.  9:pg  20] 


m'log^{N)  +  2NAr'  +  8Nlog,  (*V')  +  12N  (15) 


FSM  Complexity  [Ref.  9:pg  16] 

+  2Nlog2(N)  (16) 

The  complexity  of  four  example  cases  with  varying  N  and  N' 
are  estimated  below  using  equation  (14) ,  (15)  and  (16) . 


N 

N' 

FAM 

SSCA 

FSM 

1024 

64 

5.3*10* 

3.6*10* 

10® 

2048 

128 

2.0*10® 

1.4*10® 

4.2*10® 

32,768 

512 

1.5*10* 

1.1*10* 

10’ 

1,048,576 

8192 

8.0*10'° 

7.0*10'° 

10'^ 

It  is 

evident 

that  FAM  and 

SSCA  require 

substantially 

fewer  multiplies  than  FSM.  Particularly  with  larger  values  of 
N,  FAM  and  SSCA  require  NIogjN  while  FSM  requires  floating 
point  multiplies. 


VI.  CONCLUSIONS 


A.  SUMMARY 

The  purpose  of  this  thesis  was  to  implement  the  FAM  and 
SSCA  for  use  in  cyclic  spectral  analysis  research.  Both 
algorithms  were  successfully  implemented  in  a  manner 
compatible  with  the  SSPI  analysis  package  to  facilitate  future 
work.  By  allowing  the  user  a  choice  of  two  output  file 
formats,  results  may  easily  be  used  in  either  MATLAB  [Ref.  6] 
or  SSPI  [Ref.  4]  for  further  signal  processing  applications. 

It  has  been  shown  that  the  FAM  and  SSCA  algorithms 
generate  results  consistent  with  FSM  but  require  substantially 
fewer  floating  point  multiplies  than  FSM.  This  is  especially 
true  as  the  number  of  sample  data  points  used  increases. 

B .  RECOMMENDATIONS 

In  the  course  of  this  thesis  it  has  become  apparent  that 
the  FAM  and  SSCA  algorithms  consist  of  inherently  parallel 
operations.  Roberts  and  Loomis  explore  these  possibilities  in 
Reference  10.  The  high-level  sequential  nature  of  these 
algorithms  also  lend  them  to  pipelining  architectures.  By 
combining  a  parallel  or  multi-processor  approach  with 
pipelining  a  large  improvement  in  performance  should  be 
obtainable.  This  would  be  espectially  true  in  applications 


requiring  large  sets  of  signal  data  samples  or  repeated, 
time- sequenced  cyclic  spectral  plane  estimations. 

Two  areas  which  could  benefit  from  further  work  are 
automatic  signal  identification  and  signal  parameterization 
algorithms.  The  results  from  cyclic  spectral  analysis 
algorithms  such  as  FAM  and  SSCA  may  be  input  into  Linear 
Associator  or  Mapping  neurocomputer  networks  as  described  in 
Reference  ll.  Intensity  transformations  using  digital  image 
processing  techniques  from  Reference  12  also  appear  to  have 
promise.  Both  classes  of  techniques  could  have  utility  in 
automatic  signal  identification  and  parameterization.  They 
also  have  the  added  benefit  that  they  lend  themselves  to 
parallel  and  pipelined  computing  architectures. 
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APPENDIX  A.  FAM  PROGRAM  USE 


Correct  coiranands  are: 
for  the  cross- fam 

fam  inputfilel  inputfile2  outputfile  n  nprime  1  Oflag 
fam  iroutfilel  inputfile2  outputfile  n  nprime  1  Oflag  red 

or  for  the  auto -fam 

fam  inputfilel  outputfile  n  nprime  1  Oflag 
fam  inputfilel  outputfile  n  nprime  1  Oflag  red 

where:  inputfilel  is  the  file  containing  the  x  signal  samples 

inputfile2  is  the  file  containing  the  y  signal  samples 
outputfile  is  where  the  spectrum  values  will  be  placed 
n  is  the  number  of  samples  to  use  from  the  inputfiles 
and  it  must  be  a  power  of  2 
nprime  is  the  group  size  of  the  input  datasets  and  it 
must  be  a  power  of  2 

1  is  the  starting  offset  of  subsequent  datasets  and 
it  must  be  a  power  of  2 

Oflag  indicates  the  output  filetype,  ascii  or  binary 
-asc  for  an  ascii,  MATLAB  compatible  file, 
full  cyclic  spectral  plane 
-plo  for  a  plot_sxaf,  SSPI  compatible  file, 
half  cyclic  spectral  plane 
red  reduces  the  amount  of  data  output 

note:  the  first  two  entries  in  every  input  file  are  expected 

to  be  the  datatype  and  n.  datatype  =  l  for  real 
values,  2  for  complex  values.  n  is  the  number  of 
samples  contained  in  the  file,  one  sample  per  line. 

input  file  format: 

datatype  n 
sample  #1 
sample  #2 


sample  #n 
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output  file  ASCII  MATLAB  format: 


nprime  n 

output (1)  (1) 
output  (1)  (2) 


output (nprime) (n) 


output  file  SSPI  ploC_sxaf  format : 

datatype  (1  for  real,  2  for  complex) 
number_of_alphasalpha_minalpha_max 
number_of _f  reqs  f  req_minf  r eq_max 

value  of  alpha  #l 

number  of  freqs  in  #lf req_minf req_max 
spectrum  at  freq_min 


spectrum  at  freq_max 


value  of  alpha  #alpha_max 

number  of  freqs  in  #alpha_maxf req_minf req_max 
spectrum  at  freq_min 


spectrum  at  freq_max 
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APPENDIX  B.  FAM  PROGRAM  LISTING 
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»ti  iti  tft  t-tt 


:a.T. . , 


.-acre  . 


5  ir. c_aas  '  ?*"  alii^ .  r. ^ 

*  1  r.  c  -  u  a  0  ^  '  5  a  a  1  r  .  r.  ^ 
sir.a*uae  <a.atr.  .r.> 

^larlaae  "  ac-ei  carter  taesrs.SS?:  part  frt.a" 

=  ir.rluae  "  r.trr.eJ  carter  ■  thesis- SS?I  paa.- raaix .  c" 
r^ir.  (arac,  arpv) 
ir.t  argc; 
c.tar  'arav;;,- 

3CMP1.ZX  **s3,^*v3/  '^ciwr.ccnv,  **s4; 
rCN^lEX  'teT.pfft,  *yl,  'si,  "s2,  ’*y2,  'vir.dcw; 

'  X  passes  data  to, from  fft  routine 

33  holds  results  of  FFTing  s2 
y3  holds  results  of  FFTing  y2 

awnccnv  holds  downconversion  coefficients 

34  holds  correlation  multiply  results 
tempfft  passes  data  to/ from  fft  routine 
yl  holds  y  samples  from  inputfiie  2 

si  holds  X  samples  from  inputfile  1 
s2  holds  channelized  x  samples 
y2  holds  channelized  y  samples 
window  holds  Hcimm.ing  window  values 

*/ 

int  1, ], k, type, n, f ile_n, nprime, p, 1, direction,  norm; 
int  a,  b,  c,  cross, num_f , max_num_f , data_type,max_num_alf ; 
int  tempi, temp], halfp, reduce, curr_alf ; 
int  'outint; 

loat  numreai,  numimag, convfac,  mainfac, num_alf , f_min,  f_max; 
loat  alpha_min,  alpha_max, f_min_all, f_max_all; 
loat  tempreal, tempimag, bigmag, tempmag; 
loat  twopi=6. 23318530718; 
loat  'outreal; 
double  z,y; 

char  'infilel, *infile2, 'outfile; 

FILE  'ifpl, *ifp2, 'ofp; 

■'  check  for  the  correct  number  of  input  arguments 
the  correct  commands  are: 

fam  inputfilel  inputfile2  outputfile  n  nprime  1  Oflag  reduce 
or 

fam  inputfilel  inputfileE  outputfile  n  nprime  1  Oflag 
or 

fam  inputfilel  outputfile  n  nprime  1  Oflag  reduce 
or 

fam  inputfilel  outputfile  n  nprime  1  Oflag 

where:  inputfilel  is  where  the  x  signal  samples  are  expected  to  be 

inputflleE  IS  where  the  y  signal  samples  are  expected  to  be 

outputfile  IS  where  the  spectrum  values  will  be  put 

n  IS  the  number  of  samples  to  us“  from  the  inputfiles 
nprime  is  the  group  size  of  input  datasets 
1  IS  the  starting  offset  of  subsequent  datasets 
Oflag  indicates  the  output  filetype,  ascii  or  binary 

-asc  for  an  ascii,  matlab  compatible  file,  full  plane 
-plo  for  a  plot_sxaf  compatible  file,  half  plane 
reduce  is  an  option  to  reduce  the  amount  of  output 
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note:  datatype  and  n,  are  expected  to  ce  at  the  top  of  the  input  file 

datatype  =  1  for  real,  2  for  complex 
n  IS  the  number  of  samples  following 

written  by  LCDR  Nancy  J.  Carter,  USN  NPGS  Monterey  CA  01  October  1992 
to  implement  the  Frequency  Accumulation  Method  of  cyclic  spectral  analysis 
as  described  in  the  Brown/Gardner/Loomis  IEEE  paper 

20  October  1992 ....  changed  to  ascii  MATLAB  output  format 

22  October  1992 ... .added  plot_sxaf  output  format  compatible  with  SSPI 

23  November  1992 ... .added  option  to  reduce  amount  of  output 

» / 

/*  CHECK  INPUT  VALUES  */ 


/*  */ 
if  (  (argc  !=  7)&&(argc  !=  8)&&(argc  '.=  9) )  { 
printf ("fam. ...  fatal  errorXn"); 

printfC' . incorrect  number  of  calling  argument s\n" ) ; 

printf  (" . correct  formats  are:\n"); 

printf  (" .  fam  inputfilel  inputfiie2  outputfile  n  nprime  1  Of  lag  reduceXn") 

printf  (" .  fam  inputfilel  inputfiie2  outputfile  n  nprime  1  OflagXn"); 

printf  (" .  fam  inputfile  outputfile  n  nprime  1  Oflag  reduceXn"); 

printfC .  fam  inputfile  outputfile  n  nprime  1  OflagXn"); 

printfC .  whereXn"); 

printfC .  inputfilel  contains  the  x  samplesXn"); 

printfC .  inputfile2  contains  the  y  samplesXn"); 

printfC .  outputfile  will  contain  the  resultsXn")  ; 

printfC .  n  is  the  number  of  input  samples  to  useXn"); 

printfC .  nprime  is  the  group  size  of  input  datasetsXn" ) ; 

printfC .  1  is  the  offset  of  subsequent  datasetsXn"); 

printfC .  Oflag  is  the  output  file  formatXn"); 

printfC .  Oflag  =  -asc,  an  ascii  file  is  producedXn" )  ; 

printfC .  Oflag  =  -plo,  a  plot_sxaf  file  is  producedXn"); 

printfC .  reduce  is  an  option  for  reduced  amount  of  outpurXn"); 

exit  0  ; 

} 


if  (argc==7)  {  /*  this  is  an  auto-fam,  no  output  reduction  */ 

cross=0; 

infilel=argv[l] ; 
out  f ile=argv [ 2 ] ; 
n=atoi (argv [3] ) ; 
nprime=atoi (argv [4 ] ) ; 
l=atoi (argv [5 ] ) ; 
reduce=0; 

} 

if  ( (argc==8) i 4 (argv [7 ] [0] ! =' r' ) )  {  /*  this  is  a  cross-fam,  no  output  reduction  */ 

cross=l; 

inf ilel=argv [ 1 ] ; 
inf ile2=argv [2 ] ; 
outf ile=argv [ 3 ] ; 
n=atoi (argv [4] ) ; 
nprime=atoi (argv [ 5 ] ) ; 
l=atoi(argv[6] ) ; 
reduce=0; 

} 

if  ( (argc==8) 44 (argv[7] [0] ==' r' ) )  {  /*  this  is  an  auto-fam,  with  output  reduction  */ 

cross=0; 
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inf ilel=drgv [ 1 ] ; 
out  f ile=argv [2]; 
n=atoi(argv[3] )  ; 
nprime=atoi (argv [ 4  ]  )  ; 
l=atoi (argv [5] ) ; 
reduce=l; 

} 

if  (argc==9)  {  /*  this  is  a  cross-fam,  with  output  reduction  ’/ 

cross=l ; 

inf ilel=argv [ 1 ] ; 
infile2=argv [2 ] ; 
outf ile=argv [3 ] ; 
n=atoi (argv [4] ) ; 
nprime=atoi (argv [ 5 ] )  ; 
l=atoi (argv [ 6] ) ; 
reduce=l; 

} 

/*  verify  that  n,  nprime  and  1  are  powers  of  2  */ 

i=n%2; 

if  (i!=0){ 

printf (" fam. ...  fatal  errorXn"); 

printfC . calling  argument  n  is  not  a  power  of  2\n"); 

exit  ( )  ; 

} 

i=nprime%2; 
if  (i!=0){ 

printf ("fam. . . . fatal  errorXn") ; 

printfC . calling  argument  nprime  is  not  a  power  of  2\n"); 

exit  0  ; 

} 

]=1%2; 
if  <j!=0){ 

printf ("fam. ...  fatal  errorXn"); 


printfC . calling  argument  1  is  not  a  power  of  2Xn"); 

exit  ( )  ; 

} 

/*  */ 

/*  open  input  file  1,  prepare  to  get  the  x  signal  sample  data  */ 
/*  V 


ifpl  =  fopen (infilel, "r") ; 
fscanf (ifpl, "%i  %i", &type, &file_n) ; 

/*  verify  that  file_n  is  greater  than  or  equal  to  n  */ 
if  (file_n<n) { 

printf ("fam. . . . fatal  errorXn") ; 

printfC . inputfilel  does  not  contain  enough  samplesXn"); 

fclose (ifpl) ; 
exit  0  ; 

} 

/*  find  p  -  the  number  of  datasets (nprime  long)  */ 
p= ( (n-nprime) /I) ; 

/* 

/*  */ 

/*  READ  IN  DATA  SAMPLES  */ 

/*  */ 

/*  allocate  space  to  read  in  the  sample  values  */ 

/*  */ 


/ 
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3l= (COMPLEX*) calloc (n, sizeof (COMPLEX) )  ; 
if  (sl==NULL) { 

printf (" fam. ...  insufficient  space  to  allocate  sl\n"); 
exit  0  ; 

) 

if  (type==l)  { 

/  *  *  / 

/*  read  in  the  real  sample  values  */ 

/*  */ 

for  (i=0;  i  <  n;  i++) { 

fscanf (ifpl, "%e\n", snumreal) ; 
s  1  [  1  ]  .  r=nuinreal  ; 
si [i] . 1=0 .0; 

} 

else  { 

/*  */ 

/*  read  in  the  complex  sample  values  */ 

/  *  *  / 

for  (i=0;  i  <  n;  i++)  { 

fscanf (ifpl, "%e  %e\n", snumreal, &numimag) ; 

si [i] . r=numreal; 
si [i] . i=numimag; 

} 

} 

/*  */ 

/*  close  the  input  file  #1  */ 

/*  */ 

fclose (ifpl) ; 

/*  */ 

/*  if  this  is  a  cross-fcim  go  get  y  samples  */ 

/*  */ 

if  (cross==l)  { 

/*  */ 

/*  open  inputfile2,  prepare  to  get  the  y  signal  sample  data  */ 

/*  */ 

ifp2  =  fopen (infile2, "r") ; 
fscanf (ifp2, "%i  %i", itype, 4file_n) ; 

/*  verify  that  file_n  is  greater  than  or  equal  to  n  */ 
if  (file_n<n) { 

printf ("fam. ...  fatal  error\n"); 

printf (" . input filel  does  not  contain  enough  samples \n" ) ; 

fclose (ifpl) ; 
exit  0  ; 

} 

/*  */ 

/♦  READ  IN  Y  DATA  SAMPLES  */ 

/*  */ 

/*  allocate  space  to  hold  the  sample  values  */ 

/*  */ 

yl= (COMPLEX*) calloc (n, sizeof (COMPLEX) ) ; 
if  (yl=-NULL) { 

printf  ("fcim.  ...  insufficient  space  to  allocate  yl\n"); 
exit  0  ; 

) 

/*  */ 
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/  read  ic  the  complex  sample  values  ’ 

/•  *  »  '■ 
for  (1=0;  1  <  n;  i++)  f 

fscanf (ifp2, " oe  %e\n", inumreal, inumimag) ; 
yl [ 1 1 . r=numreal; 
y  1  [  1  ]  .  i=n\imimag; 

} 

/*  */ 

/*  close  the  inputfile2  */ 

/♦  */ 

f close (ifp2) ; 

}  /*  end  if  cross=l  */ 

/*  ♦/ 

/*  FORM  DATASETS  */ 

/*  */ 

/*  allocate  space  to  hold  the  p-by-nprime  datasets  */ 

/♦  */ 

s2= (COMPLEX**) calloc (p,  sizeof (COMPLEX*) ) ; 
if  (S2==NULL) { 

printf ("fam. .. .insufficient  space  to  allocate  s2\n"); 
exit  () ; 

} 

for  (1=0;  i<p;  i++) { 

s2 [i]= (COMPLEX*) calloc (nprime, sizeof (COMPLEX) ) ; 
if  (s2 (i]==NULL) { 

printf ("fam. ...  insufficient  space  to  allocate  s2\n”); 
exit  0 ; 

} 

/*  */ 

/*  form  p  datasets  of  nprime  samples  from  the  original  sample  stream 
/*  */ 

for  (i=0;  i<p;  i++) ( 

for  (j=0;  3<nprime;  :++) ( 
k=i*l+j; 
s2[i]  []]=sl[k]; 

} 

} 

/*  */ 

/*  if  this  is  a  cross-fcim  form  y  datasets  */ 

/*  */ 

if  (cross==l)  ( 

/*  */ 

/*  */ 

/*  allocate  space  to  hold  the  p-by-nprime  datasets  */ 

/*  */ 

y2= (COMPLEX**) calloc (p, sizeof (COMPLEX*) ) ; 
if  (y2==NULL) ( 

printf (" fam. ...  insufficient  space  to  allocate  y2\n"); 
exit  () ; 

} 

for  (i=0;  i<p;  i++) { 

y2 [i]  =  (COMPLEX*) calloc (nprime, sizeof (COMPLEX) ) ; 
if  (y2 [i]==NULL) { 

printf ("fam. ...  insufficient  space  to  allocate  y2\n"); 
exit  ( ) ; 
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form  p  dacasers  of  nprime  samples  from  che  original  sample  scream 

/  -k  k  / 

for  (1=0;  i<p;  i++) { 

for  (1=0;  :<nprime;  i++) { 
k=i*l+i ; 
y2[i] [i]=yl[k]; 

} 

} 

}  /*  end  if  cross=l  " / 

i  *  */ 

/*  APPLY  WINDOW  TO  DATASETS  */ 

/*  *  / 

/*  allocate  space  to  hold  the  window  multiplicand  values  */ 

/*  */ 

window= (COMPLEX*) calloc (nprime, sizeof (COMPLEX) ) ; 
if  {window==NULL) { 

printf ("fam. ...  insufficient  space  to  allocate  sl\n"); 
exit  ( )  ; 

) 

/*  »/ 

/*  calculate  the  window  values  for  the  "nprime"  sample  wide  window  */ 

/*  *•/ 

for  (i=0;  i<npriine;  i++)  { 
y= (twopi*i) / (nprime-1) ; 

2=cos (y) ; 

window[i].r=  0.54  -  (0.46*2); 
window [i] .1=  window [i] .r; 

} 

/*  */ 

/*  apply  the  window  to  all  rows  of  s2  */ 

/*  */ 

for  (1=0;  i<p;  i++) { 

for  (3=0;  3<nprime;  3++) { 

s2 [i]  [3] .r=s2  [i] [ j] .r*window[3] .r; 
s2 [1] [3] .i=s2 [1] [3] .i*window[3] .1; 

} 

} 

/*  */ 

/*  if  this  is  a  cross-fam  apply  window  to  y2  */ 

/*  */ 

if  (cross==l)  { 

/*  */ 
for  (i=0;  i<p;  i++) { 

for  (3=0;  3<nprime;  3++) { 

y2 [i] [3] .r=y2 [i] [ j] .r*window[3] .r; 
y2[i] [3] .i=y2[i] [3] .i*windowf3] .1; 

} 

} 

}  /*  end  if  cross=l  */ 

/*  */ 

/*  FFT  EACH  DATASET  ROW  */ 

/*  */ 

/*  allocate  space  to  hold  rows  of  data  for  passing  to  the  FFT  routine  */ 
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;<=  (COMPLEX*)  calloc  (nprime,  sizeof  (COMPLEX)  )  ; 
if  (x==NULL) i 

prir.tf  ("fam.  ...  insufficient  space  to  allocate  x\n"); 
exit  ()  ; 

} 

/■’  ’/ 

/*  allocate  space  to  hold  rows  of  results  from  the  EFT  routine 

/  *  *  / 

s3= (COMPLEX**) calloc (p,  sizeof (COMPLEX*) ) ; 
if  (s3==NULL) ( 

printf ("fam. ...  insufficient  space  to  allocate  s3\n"); 
exit  0  ; 

t 

for  (1=0;  i<p;  i++) { 

s3 [i] = (COMPLEX*) calloc (nprime, sizeof (COMPLEX) ) ; 
if  (s3 [i]==NULL) t 

printf (" fam. ...  insufficient  space  to  allocate  s3\n"); 
exit  ( )  ; 

) 

} 

/*  */ 

/*  take  an  FFT  of  each  row  of  s2  */ 

/*  */ 

/*  set  values  in  arguments  sent  to  fft  */ 

direct ion=l; 

norm«l; 

for  (i=0;  i<p;  i++) { 

/♦  copy  rcw  of  s2  into  complex  array  */ 

for  (j=0;  j<rprime;  j++) { 
x[ j]=s2 [i] [ j] ; 

} 

/*  go  get  FFT  performed  */ 

fft (x, nprime, direction, norm) ; 

/*  copy  X  va  ues  into  a  row  of  s3  */ 

for  (j=0;  ]<i. prime;  ]++)  { 
s3[i] [:]=x:  j] ; 

} 

} 

/*  */ 

/*  if  this  IS  a  cross-fam  take  an  FFT  of  each  row  of  y2  */ 

/*  */ 

if  (cross==l)  { 

/*  */ 
/*  */ 

/*  allocate  space  to  hold  rows  of  results  from  the  FFT  routine 

/♦  */ 

y3= (COMPLEX**) calloc (p,  sizeof (COMPLEX*) )  ; 
if  (y3==NULL) { 

printf ("fam. ...  insufficient  space  to  allocate  y3\n"); 
exit  0  ; 

} 

for  (i=0;  i<p;  i++) ( 

y3 [i]  =  (COMPLEX*) calloc (nprime, sizeof (COMPLEX) )  ; 
if  (y3 [i]==NULL) ( 

printf ("fam. ...  insufficient  space  to  allocate  y3\n"); 
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exit  ( )  ; 


for  (1=0;  i<p;  i-i-+)  { 

/*  copy  row  of  y2  into  complex  array  */ 

for  (]=0;  ]<nprime;  ]++) { 

X [ : ] =y2 [ 1 1  [  ]  ] ; 

} 

/*  go  get  FFT  performed  */ 

f ft (x, nprime,  direction,  norm)  ; 

/*  copy  X  values  into  a  row  of  y3  */ 

for  (]=0;  ]<nprime;  3++) { 
y3 [i] [:]=x[ j] ; 

}  /*  end  if  cross=l  */ 

/*  */ 

/*  */ 

/*  DOWNCONVERT  EACH  ROW  »/ 

/*  */ 

/*  allocate  space  to  hold  the  downconversion  multiplicands  */ 

/*  */ 

dwnconv= (COMPLEX**) calloc (p, sizeof (COMPLEX*) )  ; 
if  (dwnconv==NULL) { 

printf ("fam. ...  insufficient  space  to  allocate  dwnconv\n"); 
exit  0 ; 

} 

for  (i=0;  i<p;  i++) { 

dwnconv [i] = (COMPLEX*) calloc (nprime, sizeof (COMPLEX) ) ; 
if  (dwnconv (i ]  ==mJLL)  { 

printf ("fam. ...  insufficient  space  to  allocate  dwnconv\n"); 
exit  0  ; 

} 

} 

/*  */ 

/*  downconvert  each  of  the  transform  sequences  */ 

/*  */ 

/*  calculate  the  down  conversion  multipliers  */ 

/*  */ 

mainfac=twopi*l /nprime; 
for  (i=0;  i<p;  i++) ( 

for  (3=0;  3<nprime;  3++) { 
convfac=i*3*mainfac; 
dwnconv [ 1 ] [ 3 ] . r=cos (convfac)  ; 
dwnconv [i][3].i=  (-1.0) * sin (convfac) ; 

) 

} 

/*  multiply  downconversion  factor  against  each  frequency  value  */ 
for  (i=0;  i<p;  i++) { 
for  (3=0;  3<nprime;  3++) { 

numreal=s3 [i] [j] .r*dwnconv[i] [jj.r  -  s3[i][j] . i*dwnconv [1 ] [3] .1; 
s3 [i]  [ 3 ] . i=s3  [i]  [ 3 ] .r*dwnconv [i]  [3].!  +  s3[i][3] . i*dwnconv [i]  [ 3 ] . 
s3[i] [j] .r=numreal; 

) 

} 

if  (cross==l)  { 
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*  ’  / 

/  if  this  13  a  cross-fam  downconvert  y3  * 

/*  */ 

/*  multiply  downconversion  factor  against  each  frequency  value  '/ 
for  (1=0;  i<p;  i++) { 

for  (3=0;  3<nprime;  ]++) { 

numreal=y3 [i] [ ] ] . r’dwnconv [i] [j].r  -  y3[i][j] . I'dwnconv [i ] [ 3 ] . 1; 
y3  [i]  [  3  ]  .  i=y3  [1  ]  [  3  ]  .  r*dwnconv  [i][3].i  ■*-  y3[i][3].  I'dwnconv  [1  ]  [  -  ]  .  r; 
y3[i] [3] . r=numreal; 

} 

}  /*  end  if  cross=l  */ 

/  *  *  / 

/»  MULTIPLY  COLUMNS  */ 

/*  */ 

/*  allocate  space  to  hold  the  correlation  mult^-piy  (column  multiply)  results  */ 
/*  */ 

a=nprime*p/2; 

s4= (COMPLEX**) calloc (a, sizeof (COMPLEX*) ) ; 
if  (s4==NULL) { 

printf ("fam.  ...  insufficient  space  to  allocate  s4\n"); 
exit  0  ; 

} 

for  (i=0;  i<a;  i++) { 

s4 [i] = (COMPLEX*) calloc (nprime, sizeof (COMPLEX) ) ; 

If  (s4[i]==NULL) { 

printf ("fam. ...  insufficient  space  to  allocate  s4\n”); 
exit  0  ; 

} 

} 

/*  */ 

/*  allocate  space  to  hold  columns  of  s4  for  passing  to  the  FFT  routine  */ 

/*  */ 

teinpfft=  (COMPLEX*)  calloc  (p,  sizeof  (COMPLEX) )  ; 
if  (tempfft==NULL) ( 

printf ("fam. .. .insufficient  space  to  allocate  tempfft  %i  by  l\n",p); 
exit  0  ; 

1 

/*  set  values  in  arguments  sent  to  fft  */ 

direction=l; 

norm=l; 

if  (cross==0)  {  /*  auto-fam  */ 

/*  */ 

/*  multiply  columns  of  s3  */ 

/*  */ 

for  (a=nprime-l;  a>  -1;  a — )  { 
for  (b=0;  b<nprime;  b++)  { 

/*  multiply  the  columns  of  s3  with  other  columns  of  s3  */ 
for  (3=0;  3<p;  3++) { 

tempfft [ j] .r= (s3 [3] [a] .r*s3 [3] [b] .r)+(s3 [ j] [a] .i*s3 [3] [b] . 1) ; 
tempfft [3] .i=(s3[3] [a] .i*s3[jl [b] .r)-(s3(3] [a] .r*s3[j] [b] .1); 

) 

/*  */ 

/*  FFT  EACH  COLUMN  */ 

/*  */ 

/*  go  get  FFT  performed  */ 


70 
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f ft (tempf ft, p, direction, norm)  ; 

/*  copy  tempfft  result  values  back  into  a  row  of  s4  ’'/ 

for  (2=(p/4);  ]< ( (3*p/4) -1)  ;  ]++)  { 
c=  (  (nprime-a-1 )  *p/2 )  -r^-  (p/  4)  ; 
s4[c] [b] =tempf ft [ ] ] ; 

} 

}  / *  end  for  b=0 ...  */ 

}  / *  end  for  a=0 ...  * / 

}  /*  end  auto-fam  */ 

else  {  /*  cross-fam  */ 

/  *  ♦/ 

/*  multiply  columns  of  53  with  y3  */ 

/*  */ 

for  (a=nprime-l;  a>  -1;  a — )  { 
for  (b=0;  b<nprime;  b++)  { 

/*  multiply  the  columns  of  s3  with  other  columns  of  y3  */ 
for  (j=0;  :<p;  ]++) { 

tempfft [:] .r=  (s3  []  ]  [a] .r’y3 []]  [b] . r)  +  (s3 [ ] ] [a] . i’y3  [] ]  [b] . i) ; 
tempfft []] .i=(s3[:]  [a]  .i*y3 []] [b] . r) - (s3 [ ] ]  [a]  .r’y3 [ j]  [b] .i) ; 

} 

/*  */ 

/*  FFT  EACH  COLUMN  */ 

/*  */ 

/*  go  get  FFT  performed  */ 

fft (tempfft, p,  direction,  norm) ; 

/*  copy  tempfft  result  values  back  into  a  row  of  s4  */ 

for  (j=(p/4);  j< ( (3*p/4) -1) ;  ]++)  { 
c= ( (nprime-a-1) *p/2) +j- (p/4) ; 
s4[c] [b] =tempf ft [ j ] ; 

} 

}  /*  end  for  b=0 ...  */ 

)  /*  end  for  a=0...  */ 

}  /*  end  cross-fam  ♦/ 

/*  */ 

/*■  OUTPUT  */ 

I*  *  / 

halfp=p/2; 

if  (reduce==0)  {  /*  do  not  reduce  the  amount  of  output  */ 

if  ( ( (argc==8) && (argv[7 ] [l ] ==' a' ) ) I i ( (argc==7) && (argv [6] [1 ] ==’ a'  ) ) )  ( 
/*  make  an  ascii  output  file  of  the  whole  spectral  plane  */ 

/*  open  the  output  file  and  place  header  information  in  it  */ 
ofp  =  fopen (outf ile, "w" ) ; 
max_num_alf=nprime*p/2 ; 
max_niam_f =npr  ime ; 

fprintf (ofp, "%i  %i\n",max_num_aif,max_num_f ) ; 
for  (i=0;  i<max_num_alf ;  i++) { 
for  (j=0;  ]<nprime;  ]++) { 

fprintf  (ofp,  "%e  %e\n",  s4  [i]  [;]]  .r,  s4  [i]  [3]  .  1)  ; 

}  /*  end  for  ]= . . .  */ 

}  /*  end  for  i=. . .  * / 

fclose (ofp) ; 

}  /•  end  if  ascii...  */ 

else  {  /*  make  a  plot_sxaf  no  reduced  output  file  */ 
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*  open  the  output  file  and  place  header  informaticn  in  it  ■'/ 
ofp  =  fopen (outf lie, "w" ) ; 

/*  place  header  information  in  the  file  */ 
data_type=2;  /*  real_or_complex  ' 
max_num_alf=nprime*p/2 ;  /’  num_alpha  */ 

alpha_min=0 . 0 ;  /*  alpha_min  */ 

alpha_max=l .  0;  /■*  alpha_max  ’/ 

max_num_f=nprime;  /*  Tnax_num_f  *•' 
f_min_all=  -0.5;  /'  f_min_all  ’/ 

f_max_all=  0.5;  /*  f_max_all  */ 

fprintf (ofp, "%i  %i  %e  %e  %i  %e  %e\n”, data_type, max_num_alf ,  alpha_min, 
alpha_max,  max_niim_f ,  f_min_all,  f_max_all)  ; 
for  (i=0;  i<nprime;  i++) {  /*  group  of  lines  */ 

for  (j=0;  j<(p/2);  ]++) {  /*  line  in  the  group  */ 

/*  place  alpha  line  header  information  in  file  */ 
curr_alf= (i*p/2) +];  /*  alpha  number  */ 

num_f=nprime-i;  /*  num_f  */ 

f_min=  -0 . 5+ (0 . 5*i/nprime) ;  /*  f_min  */ 
f_max=  0 . 5- (0 . 5*i/nprime) ;  /*  f_max  *! 
fprintf (ofp, " %i  %i  %e  %e\n", curr_alf , num_f , f_min,  f_max)  ; 
for  (k=i;  k<nprime;  k++) {  /*  points  on  the  line  */ 

tempi= (p/2* (nprime-1) ) - ( (k-i) *p/2) + j; 
tempj=k; 

numreal=s4 [tempi] [temp]] .r; 
numimag=s4 [tempi] [tempj] .i; 
fprintf (ofp, "%e  %e\n", numreal, numimag) ; 

}  /*  end  for  k=. . .  */ 

}  /*  end  for  j=. . .  */ 

}  /*  end  for  i=. . .  */ 

/*  close  the  output  file  */ 
f close (ofp) ; 

}  /*  end  else  ...  */ 

}  /*  end  if  reduce=0  ...  */ 

else  {  /*  reduce  the  amount  of  output  */ 

if  ( ( (argc==9)&& (argv[7]  [l]=='a' ) ) I  I ( (argc==8)&& (argv[6] [l]=='a'  ) ) )  { 

/*  make  an  ascii  output  file  of  the  whole  spectral  plane  */ 

/*  open  the  output  file  and  place  header  information  in  it  */ 
ofp  =  fopen (outfile, "w") ; 
max_num_alf =nprime ; 
max_num_f =npr ime ; 

fprintf (ofp, "%i  %i\n", max_num_alf , max_num_f )  ; 
for  (i=0;  i<max_num_alf ;  i++) { 
for  (j=0;  j<nprime;  j++) { 
numreal=s4 [i*halfp] [j] .r; 
numimag=s4 [i*halfp]  [j]  .i; 

bigmag=sqrt (numreal*numreal  +  numimag*numimag) ; 
for  (k=l;  k<halfp;  k++)  {  /*  look  for  largest  value  */ 

tempreal=s4 [i*halfp+k] [ j] .r; 
tempimag=s4 [i*halfp+k] [j] .i; 

tempmag=3qrt (tempreal*tempreal  +  tempimag*tempimag) ; 
if  (tempmag>bigmag)  (  /»  found  a  larger  value  */ 

b igmag- t empmag ; 
numreal=tempreal; 
numimag=t  emp imag ; 
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}  /♦  end  if  tempmag>bigmag.  .  .  " 

}  / *  end  if  k.=l  ...  ’ / 

fprincf  (ofp,  "le  %e \n" ,  numreai,  .numimag) ; 

}  / *  end  for  ]= . . .  * / 

}  / *  end  for  i=  .  .  .  ' / 

f close (ofp) ; 

\  /*  end  if  ascii...  */ 

else  {  /*  make  a  plot_sxaf  reduced  output  file 

/*  open  the  output  file  and  place  header  information  i.n  it  */ 
ofp  =  fopen (outf ile, "w" ) ; 

/*  place  header  information  in  the  file  */ 
data_type=2;  /*  real_or_complex  */ 
max_num_alf=nprime;  /*  num_aipha  */ 
alpha_min=0 . 0;  /*  alpha_min  ’/ 

alpha_max=l . 0;  /*  alpha_max  */ 

max_num_f=nprime;  /*  max_num_f  */ 
f_min_all=  -0.5;  /*  f_min_all  */ 

f_max_all=  0.5;  /*  f_max_all  »/ 

fprintf (ofp, "%i  %i  %e  %e  %i  %e  %e\n" , data_type, max_num_alf ,  alpha_min, 
alpha_max, max_num_f , f_min_ail, f_max_all) ; 
for  (i=0;  i<nprime;  i++) {  /*  group  of  lines  */ 

/*  place  alpha  line  header  information  in  file  */ 
curr_alf=i;  /*  alpha  */ 
num_f=nprime-i;  /*  num_f  */ 
f_min=  -0 . 5+ (0 . 5*i/nprime) ;  /*  f_min  */ 

f_max=  0 . 5- (0 . 5*i/nprime) ;  /*  f_max  */ 

fprintf (ofp, "%i  %i  %e  %e\n",curr_alf,n;un_f, f_min,  f_max)  ; 
for  (k=i;  k<nprime;  k++) {  /*  points  on  the  line  */ 

tempi® (halfp* (nprime-1) ) -  ( (k-i) *halfp) ; 
tempj=k; 

numreal=s4 [tempi] [tempj] .r; 
numimag=s4 [tempi] [tempj] .i; 

bigmag=sqrt  (numreal*n'umreal  +  numimag*numiinag) ; 
for  (j=l;  j<halfp;  ]++) [  /*  look  for  the  largest  value  */ 

tempreal=s4 [tempi+ j ] [temp^] .r; 
tempimag=s4 [tempi+ j ] [temp^] .i; 

tempmag=sqrt (tempreal*tempreal  +  tempimag*tempimag) ; 
if  (tempmag>bigmag) {  /*  found  a  larger  value  */ 

b igmag=t  empmag ; 
numreal=terapreal; 
numimag=t  emp imag ; 

}  /♦  end  if  tempmag>  ....  */ 

}  / ♦  end  for  j= . . .  */ 

fprintf (ofp, "%e  %e\n", numreal, numimag) ; 

}  /*  end  for  k=. . .  * / 

}  /*  end  for  i=. . .  */ 

/*  close  the  output  file  */ 
fclose (ofp) ; 

}  /*  end  else  binary...  ♦/ 

}  /*  end  reduce  the  amount  of  output  */ 

]  /*  end  of  mam  */ 
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Correct  commands  are: 
for  the  cross-ssca 

ssca  inputfilel  inputfile2  outputfile  n  nprime  1  Oflag 
ssca  inputfilel  inputfile2  outputfile  n  nprime  1  Oflag  red 

or  for  the  auto-ssca 

ssca  inputfilel  outputfile  n  nprime  1  Oflag 
ssca  inputfilel  outputfile  n  nprime  1  Oflag  red 

where:  inputfilel  is  the  file  containing  the  x  signal  samples 

inputfile2  is  the  file  containing  the  y  signal  samples 
outputfile  is  where  the  spectrum  values  will  be  placed 
n  is  the  number  of  samples  to  use  from  the  inputfiles 
and  it  must  be  a  power  of  2 
nprime  is  the  group  size  of  the  input  datasets  and  it 
must  be  a  power  of  2 

1  is  the  starting  offset  of  subsequent  datasets  and 
it  must  be  a  power  of  2 

Oflag  indicates  the  output  filetype,  ascii  or  binary 
-asc  for  an  ascii,  MATLAB  compatible  file, 
full  cyclic  spectral  plane 
-plo  for  a  plot_sxaf,  SSPI  compatible  file, 
half  cyclic  spectral  plane 
red  reduces  the  amount  of  data  output 

note  1:  the  first  two  entries  in  every  input  file  are 

expected  to  be  the  datatype  and  n.  datatype  =  1  for  real 
values,  2  for  complex  values.  n  is  the  number  of  samples 
contained  in  the  file,  one  sample  per  line. 

note  2:  error  may  occur  using  the  reduce  option  if 

(n- nprime) /nprime  is  not  an  integer. 

input  file  format: 

datatype  n 
sample  #1 
sample  #2 


sample  #n 
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output  ASCII  MATLAB  format: 

nprime  n 

output (1) (1) 
output (1) (2) 


output (nprime) (n) 


output  SSPI  plot_sxaf  format: 

datatype  (1  for  real,  2  for  complex) 
number  of  alphasalphaminalphamax 
number  of  f reqsf reqminf reqmax 

value  of  alpha  #1 

number  of  freqs  in  #lf reqmin_in_#lf reqmax_in_#l 
output  for  f reqmin_in_#l 


output  for  f reqmax_in_#l 


value  of  alpha  #alphamax 

number  of  freqs  in  #alphainaxf  reqminf  reqmax 
output  for  f reqmin_in_#alphamax 


output  for  f reqmax_in_#alphamax 
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^include  <stdlib.h> 

4 include  <stdio.h> 

#include  <raath.h> 

4 include  " /home! /carter/ thesis/ SSPI/pam/f ft . c" 

4 include  " /home3/carter/thesis/SSPI/pam/ radix . c" 
mam  (argc,  argv) 
int  argc; 
char  *argv[]; 

{ 

COMPLEX  *x, **s3, **dwnconv; 

COMPLEX  ’tempfft,  *sl,  *yl,  ■'*s2,  ‘window; 
int  1,  j,  k, type, n, nprime, p, i, r, direction, norm; 
int  a,  cross,  f ile_n,  m,  reduce, si2es3,  redindex; 
int  tempi, temp j, tempk, tempi; 
int  ‘outint; 

float  numreal, numimag, convfac,  mainfac; 
float  bigmag, tempmag, tempreal, tempimag; 
float  twopi=6. 28318530718; 
float  ‘outreal; 
double  z,y; 

char  *infilel, *infile2, ‘outfile; 

FILE  *ifpl, •ifp2, *ofp; 

/*  X  passes  data  to/from  fft  routine 

s3  holds  results  of  operations  after  first  fft 
si  holds  the  x  sample  values  from  input  file  1 

yl  holds  the  y  sample  values  from  inputfile  2 

s2  holds  the  channelized  x  input  data  from  si 

*  / 

/* 

check  for  the  correct  number  of  input  argiunents 
the  correct  commands  are  : 

ssca  inputfile  outputfile  n  nprime  1  Oflag 
ssca  inputfile  outputfile  n  nprime  1  Oflag  reduce 
ssca  inputfilel  inputfile2  outputfile  n  nprime  1  Oflag 

ssca  inputfilel  inputfile2  outputfile  n  nprime  1  Oflag  reduce 

where:  inputfilel  is  where  the  x  signal  samples  are  expected  to  be 

inputfile2  is  where  the  y  signal  samples  are  expected  to  be 
outputfile  IS  where  the  spectrum  values  will  be  put 
n  IS  the  number  of  input  samples  to  use 
nprime  is  the  group  size  of  input  datasets 
1  is  the  starting  offset  of  subsequent  datasets 
Oflag  indicates  the  output  file  type 

-asc  for  an  ascii,  matlab  compatible  file 
-plo  for  an  ascii,  plot_sxaf  compatible  file 
reduce  is  an  option  to  generate  a  reduced  amount  of  output 

note  1:  type  and  file_n  are  expected  to  be  on  the  first  line  of  the  inputfile/s 
type  =  1  for  real,  2  for  complex, 

file_n  =  niimber  of  samples  following  the  first  line 

note  2:  error  may  occur  if  (n-nprime)  is  not  evenly  divisible  by  nprime 
when  using  the  reduce  option 


maintenance  history 

. . . .written  by  LCDR  Nancy  J.  Carter,  USN  NPGS  Monterey  CA  01  October  1992 


?7 


Dec  9  12:22  1992  strip/ssca.c  Page  2 


....15  Oct  92 ....  copied  from  current  version  of  fam.c  to  alter  to  ssca  aigorit.hm 
. . .modified  to  perform  ssca  per  Brown/Gardner/ Loomis  paper 
....20  Oct  92.... made  outputfile  ascii  format  MATLAB  compatible 
....22  Oct  92.... made  an  output  format  3SPI  plot_sxaf  compatible 
....24  Nov  92.... added  option  for  a  reduced  amount  of  output 
V 

/  *  ’  / 

/*  CHECK  INPUT  VALUES  */ 

/  *  */ 

if  ( (argc  !=  7)&&(argc  !=  8)&&(argc  !=  9) ) { 
printf  (’’ssca.  .  .  fatal  error\n”)  ; 

printf(" . incorrect  number  of  calling  arguments \n" ) ; 

printf  (" . correct  formats  are:\n"); 

printf  (" .  ssca  inputfile  outputfile  n  nprime  1  Oflag\n"); 

printf  (" .  ssca  inputfile  outputfile  n  nprime  1  Oflag  reduce\n"); 

printf  (" .  ssca  inputfilel  inputfile2  n  outputfile  nprime  1  OflagXn"); 

printf (" .  ssca  inputfilel  inputfile2  n  outputfile  nprime  1  Oflag  reduceVn") 

printf  (" .  whereVn"); 

printf (" .  inputfilel  contains  the  x  signal  samplesVn" ) ; 

printf (" .  inputfile2  contains  the  y  signal  samplesXn”); 

printf (" .  outputfile  will  contain  the  resultsXn" ) ; 

printf  (" .  n  is  the  number  of  input  samples  to  use\n"); 

printf  (" .  nprime  is  the  group  size  of  input  datasetsXn" )  ; 

printf  (" .  1  is  the  offset  of  subseqiient  datasetsXn”); 

printf  (" .  Oflag  is  the  output  file  formatXn"); 

printf  (" .  Oflag  =  -asc,  a  MATLAB  file  is  producedXn" )  ; 

printf (" .  Oflag  =  -plo,  a  plot_sxaf  file  is  producedXn"); 

printf  (" .  reduce  is  an  option  for  reduced  amount  of  outputXn’’); 

exit  ( ) ; 

} 

if  (argc==7)  {  /♦  this  is  an  auto-ssca,  no  output  reduction  */ 

cross=0; 


inf ilel=argv [ 1 ] ; 
outfile=argv(2] ; 
n=atoi (argv[3] ) ; 
nprime=atoi (argv [ 4 ] ) ; 
l=atoi (argv [5] ) ; 
reduce=0; 

} 

if  ( (argc==8) && (argv [7 ] [0] ! =' r' ) )  {  /*  this  is  a  cross-ssca,  no  data  reduction  */ 

cross=l; 

inf ilel=argv [ 1 ] ; 
inf ile2=argv [2 ] ; 
outfile=argv[3] ; 
n=atoi (argv [4] ) ; 
nprime=atoi (argv [ 5 ] ) ; 
l=atoi (argv [6] ) ; 
reduce=0; 

} 

if  ( (argc==8) && (argv [7 ]  [0] ==' r' ) )  {  /♦  this  is  an  auto-ssca,  with  data  reduction  */ 

cross»0; 

inf ilel»argv [ 1 ] ; 
outfile“argv[2] ; 
n-atoi (argv [3] ) ; 
nprime-atoi (argv[4] )  ; 
l=atoi (argv [5] ) ; 
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reduce=l ; 


if  (argc==9)  {  /*  this  is  a  cross-ssca,  with  data  reduction  ' 

cross=l; 

inf ilel=argv [ 1 1 ; 
infile2=argv [2  I ; 
outf ile=argv [3 ] ; 
n=atoi (argv [4] ) ; 
nprime=atoi (argv [5] ) ; 
l=atoi (argv [6] ) ; 
reduce=l; 

} 

/*  verify  that  n,  nprime  and  1  are  powers  of  2  ♦/ 

i=n%2 ; 
if  (i!=0){ 

print f ( "ssca . . . fatal  error \n" ) ; 

printfC . calling  argument  n  is  not  a  power  of  2\n"); 

exit  0  ; 

} 

i=nprime%2; 

If  (i!=0)( 

printf ("ssca. . . fatal  error\n") ; 

printfC . calling  arg;iment  nprime  is  not  a  power  of  2\n"); 

exit  0  ; 

} 

j=l%2; 

If  (j!=0){ 

printf ("ssca. . . fatal  errorXn") ; 


printfC . calling  argument  1  is  not  a  power  of  2\n"); 

exit  0  ; 

} 

/*  »/ 

/*  open  input  file  1,  prepare  to  get  the  x  signal  sample  data  */ 
/*  ♦/ 


ifpl  =  fopen (infilel,  "r" ) ; 
fscanf (ifpl, "%i  %i", &type, &file_n) ; 

/*  verify  that  file_n  is  at  least  equal  to  n  »/ 
if  (file_n  <  n) { 

printf ( "ssca . . . fatal  error\n" ) ; 

printfC . inputfilel  does  not  contain  enough  samplesXn"); 

f close (ifpl) ; 
exit  () ; 

} 

/*  find  p  -  the  number  of  datasets (nprime  long)  */ 
p= ( (n-nprime) /I) ; 

/*  find  one  dimension  of  final  data  array  */ 
sizes3*p*l; 

/*  */ 

/*  READ  IN  DATA  SAMPLES  */ 

/*  */ 

/*  allocate  space  to  read  in  the  x  sample  values  */ 

/*  */ 

si- (COMPLEX* ) calloc (n,  sizeof (COMPLEX) ) ; 
if  (sl—NULL)  ( 

printf ("ssca. . .insufficient  space  to  allocate  sl\n"); 
exit  0  ; 
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if  (type==l)  {  /*  read  in  the  real  sample  values  '  ' 

for  (i=0;  i<n;  i++)  { 

fscanf (ifpl, "%e\n", snumreal) ; 
si [i] . r=numreal; 
sl[i] .1=0.0; 

} 

else  {  /*  read  in  the  complex  sample  values  */ 

for  (i=0;  i  <  n;  i++) { 

fscanf (ifpl, "%e  %e\n", snumreal, &numimag) ; 

sl[i] .r=numreal; 
sl[i] .i=numimag; 

} 

} 

/  *  */ 

/*  close  the  input  file  #1  ♦/ 

/*  */ 

fclose (ifpl) ; 

/*  */ 

if  (cross==l)  { 

/*  */ 

/*  open  input file2,  prepare  to  get  the  y  signal  sample  data  ♦/ 

/  *  •  / 

ifp2  =  fopen(infile2,  "r") ; 
fscanf (ifp2, "%i  %i", &type, &file_n) ; 

/*  verify  that  file_n  is  at  least  as  big  as  n  */ 
if  (file_n  <  n) { 
printf ("ssca. .. fatal  errorXn"); 

printfC . inputfile2  does  not  contain  enough  samplesin") 

fclose (ifp2) ; 
exit  () ; 

} 

/*  *  / 

/*  allocate  space  to  read  in  the  y  sample  values  */ 

/*  */ 

yl= (COMPLEX*) calloc (n, sizeof (COMPLEX) )  ; 
if  (yl==NULL) { 

printf ("ssca. .. insufficient  space  to  allocate  yl\n"); 
exit  0  ; 

} 

/*  ‘/ 

/*  read  in  the  complex  sample  values  */ 

/*  •/ 

for  (i=0;  i  <  n;  i++)  { 

fscanf (ifp2, "%e  %e\n", snumreal, Snumiroag) ; 

yl[i] .r=numreal; 
yl [ i ] . i=numimag; 

} 

/*  */ 

/*  close  the  inputfile2  */ 

/*  */ 

fclose (ifp2) ; 

}  /*  end  if  cross  «  1  */ 

/*  */ 

/*  FORM  DATASETS  */ 
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*  / 

- ’  allocate  space  to  hold  the  p-by-nprine  datasets  */ 

’  *  / 

s2= (COMPLEX**) calloc (p, sizeof (COMPLEX*) ) ; 

If  (s2==NULL) ( 

printf ("ssca ...  insufficient  space  to  allocate  s2\n"); 
exit  ( ) ; 


for  (1=0;  i<p;  1+^) { 

s2 [i] = (COMPLEX*) calloc (nprime, sizeof (COMPLEX) ) ; 
if  (s2 [i]==NULL) { 

printf ("ssca ...  insufficient  space  to  allocate  s2Vn"); 
exit  0  ; 

f 


/  * 

/*  form  p  datasets  of 
/* 

for  (i=0;  i<p;  i++) { 
for  (j=0;  ]<nprime; 
)c=i*l+ j ; 
s2 [i] [ j]=sl [k] ; 

} 

/*  */ 

/*  APPLY  WINDOW  TO 

/*  */ 


*  / 

nprime  samples  from  the  original  sample  stream 

*/ 


DATASETS  */ 


/♦  allocate  space  to  hold  the  window  multiplicand  values  */ 
/*  */ 


*/ 


window= (COMPLEX*) calloc (nprime, sizeof (COMPLEX) ) ; 
if  (window==NULL) { 

printf ("ssca ...  insufficient  space  to  allocate  sl\n"); 
exit  () ; 

/*  */ 

/*  calculate  the  window  values  for  the  "nprime"  sample  wide  window  */ 
/*  */ 

for  (i=0;  i<nprime;  i+-r)  { 
y= (twopi*i) / (nprime-1) ; 
z=cos (y) ; 

window[i].r=  0.54  -  (0.46*z); 
window [i] .1=  window [i] .r; 


/*  */ 

/*  apply  the  window  to  rows  of  s2  containing  data  */ 

/*  */ 

for  (i=0;  i<p;  i++)  { 

for  (j=0;  ]<nprime;  o++) { 

s2[i]  [j] .i=s2[i]  []]  .i*window[3] .i; 
s2[i]  [j] .r=s2[i] [j] .r*window[j] .r; 

} 

} 

/*  */ 

/*  FFT  EACH  DATASET  ROW  */ 

/*  */ 

/*  allocate  space  to  hold  rows  of  s2  for  passing  to  the  FFT  routine  */ 


fil 
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,  -w  »  / 

x=  (COMPLEX*)  callcc  (nprir.e,  sizeof  (COMPLEX) )  ; 
if  (x==NULL) { 

printf ( "ssca ...  insufficient  space  to  allocate  x\n"); 
exit  ( ) ; 


/  *  *  / 

/*  allocate  space  to  hold  rows  of  results  from  the  EFT  routine 

/  ★  It ! 

s 3= (COMPLEX**) calloc (nprime, sizeof (COMPLEX*) ) ; 
if  (s3==NULL) { 

printf ("ssca ...  insufficient  space  to  allocate  33\n”); 
exit  ( )  ; 

} 

for  (1=0;  i<nprime;  i++) ( 

s3 [i]= (COMPLEX*) calloc ( (p*l) , sizeof (COMPLEX) ) ; 
if  (s3 [i]==NULL) { 

printf ( "ssca ...  insufficient  space  to  allocate  s3\n"); 
exit  ( )  ; 

} 


/*  * ! 

/*  take  an  FFT  of  each  row  of  s2  */ 

/*  */ 

/*  set  values  in  arguments  sent  to  fft  */ 

direction=l; 

norm=l; 

for  (i=0;  i<p;  i++) { 

/*  copy  row  of  s2  into  complex  array  */ 

for  (j=0;  ]<nprime;  ]++) { 
x[ j]=s2 [i]  [  jl ; 

} 

/*  go  get  FFT  performed  */ 

fft (x,  nprime,  direction, norm) ; 

/*  copy  X  values  into  1th  column  of  s3  */ 

r=i*l; 

for  (]=0;  3<nprime;  :++) ( 
s3[3] [r]=x[3] ; 

} 

t 

/*  */ 

/*  DOWNCONVERT  EACH  ROW  */ 

/*  */ 

/*  allocate  space  to  hold  the  downconversion  multiplicands  */ 

/*  */ 

dwnconv= (COMPLEX**) calloc (p, sizeof (COMPLEX*) )  ; 
if  (dwnconv==NULL) { 

printf ("ssca. .. insufficient  space  to  allocate  dwnconvSn"); 
exit  0  ; 

} 

for  (i=0;  i<p;  i++) { 

dwnconv [i] « (COMPLEX*) calloc (nprime, sizeof (COMPLEX) ) ; 
if  (dwnconv [i] ==NULL) { 

printf ("ssca. . .insufficient  space  to  allocate  dwnconvXn") ; 
exit  0  ; 

} 


fiS 
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/’  lownconvert  each  of  the  transform  sequences 

,  »  *  / 

. ^  calculate  the  down  conversion  multipliers 

*/ 

mainf ac=twopi ’ 1 / npr ime ; 
for  (m=0;  m<p;  m++) { 
for  (lc=0;  jc<nprime;  k++)  { 

convfac=mainfac* (k- (nprime/2 ) +1) *m; 
dwnconv[m] [k] . r=cos (convfac) ; 
dwnconv [m] [k] . i=  (-1 .0) *sin (convfac) ; 

} 


/*  multiply  downconversion  factor  against  each  frequency  value  */ 
for  (1=0;  i<p;  i++) { 
tempi=i*l; 

for  (]=0;  ]<nprime;  ]++) ( 

numreal=s3 [ ] ] [tempi] . r 'dwnconv [i][3].r  -  s3(]] [tempi] . i 'dwnconv [ i] [ ] ] . i; 

s3  [  :  ]  [tempi  ]  .  i=s3  [  j  ]  [tempi]  .  r 'dwnconv  [i][]].i  +  s3(:j]  [tempi  ]  .  i 'dwnconv  [  i]  [  j  ]  . 

s3[3] [tempi] .r=numreal; 

} 

} 

/*  */ 

/»  REPLICATE  COLUMNS  */ 

/*  */ 

for  (i=0;  i<p;  i++)  [ 
r=i*l; 

/*  copy  column  into  1-1  adjacent  columns  */ 
for  (j=0;  ]<nprime;  j++)  { 
for  (k=l;  k<l;  k++)  { 
s3 [ j] [r+k] =s3 [ j] [r] ; 

}  /*  end  for  k=. . .  */ 

}  /*  end  for  3= . . .  * / 

}  /*  end  for  1= . . .  */ 

/*  */ 

/'  MULTIPLY  COLUMNS  */ 

/*  */ 

if  (cross==0)  {  /*  auto-ssca  '/ 

/'  ♦/ 

/'  multiply  columns  of  s3  with  conjugate  value  of  si  '/ 

/*  */ 

for  (a=0;  a<sizes3;  a++)  { 

for  (j=0;  :<nprime;  ]++) { 

s3[j] [a] .r=(s3[j] [a] .r*sl[al .r)  +  (s3[]] [a] .i*sl[a]  .1); 
s3[j] [a] .i=(s3[j] [a] .i*sl[a] .r)-(s3[]] [a] .r*sl[a] .1); 

} 

}  /'  end  for  a=0 ...  */ 

]  /*  end  auto-ssca  '/ 

else  [  /*  cross-ssca  •/ 

/*  *  / 

/*  multiply  columns  of  s3  with  conjugate  value  of  yl  */ 

/'  */ 

for  (a=0;  a<sizes3;  a++)  { 

for  (j=0;  j<nprime;  j++) { 

s3[j] [a] .r-(s3[j] [a] .r*yl[a] .r)+(s3[jl [a] .i*yl[a] .1); 
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s3[:j  [a] .i=(53[]] [ai . i*yl [a].r)-(s3[:;;a;.r'yi;ai.i); 

•  /'  ’  end  for  a=0  ...  ’  / 

.  end  crcss-ssca  " ' 

/*  FFT  EACH  ROW  ’/ 

V 

*  / 

i' "  allocate  space  to  hold  rows  of  s3  for  cassing  to  the  FFT  routine  ’/ 

/  *  *  / 

tempf ft= (COMPLEX* ) calloc (si2es3, sizeof (COMPLEX) ) ; 
if  '(tempfft==NULL)  { 

printf ("ssca. .. insufficient  space  to  allocate  tempf ft  %i  by  1  \n" , sizes3 ) ; 
exit  ( )  ; 

/*  set  values  in  arguments  sent  to  fft  */ 

direction=l; 

norm=l; 

for  (]=0;  3<nprime;  ]++) { 

/*  copy  row  of  s3  into  tempfft  */ 
for  (k=0;  tc<sizes3;  )c++)  ( 
tempfft  [k]  =s3  [  3  ]  [k.] ; 

} 

/*  go  get  FFT  performed  */ 

fft (tempfft, 3 izes3, direction, norm) ; 

/•  copy  tempfft  result  values  back  into  a  row  of  s3  »/ 

for  (k=0;  k<sizes3;  k++) { 
s3 [ j ] [k] “tempfft [k] ; 

} 

}  /*  end  for  ]=0 ...  */ 

/*  */ 

/*  OUTPUT  */ 

/*  */ 

if  (reduce==0)  {  /*  no  output  reduction  */ 

if  ( ( (argc==7)4& (argv[6]  [l]=='a' ) ) I  I ( (argc==6)&& (argv[5] [l]=='a'  ) ) )  { 

/»  make  an  ascii  output  file  full  spectral  plane  */ 

/*  open  the  output  file  and  place  header  information  in  it  */ 

ofp  =  fopen (outf lie, "w" ) ; 

fprintf (ofp, "%i  %i\n", nprime, sizes3)  ; 

/*  write  out  d^ta  values  to  file  for  plotting  */ 
for  (1=0;  i<nprime;  i++) { 
for  (j=(;  j<sizes3;  j++) { 

fprintf (ofp, "%e  %e\n",s3[i] [3] .r,s3[il [3] .1); 

}  /*  for  3=0 ...  */ 

)  / *  for  i=0 ...  */ 

/*  close  the  output  file  */ 
fclose (ofp) ; 

}  /♦  if  ascii ...  */ 

else  {  /*  make  a  plot_sxaf  output  file  half  spectral  plane  */ 

/*  open  the  output  file  and  place  header  information  in  it  */ 
ofp  =  fopen (outf ile, "w" ) ; 
a-0;  /♦  alpha  =  0  */ 

fprintf (ofp, "%i  %i  %e  %e  %i  %e  %e\n", 2, sizes3, 0 . 0, 1 . 0, nprime, -0 . 5,  0 . 5) ; 
for  (i-0;  i<nprime;  i++)  {  /♦  group  of  lines  */ 

for  (j-0;  3< (sizes3/nprime) ;  3++)  {  /*  line  in  the  group  */ 

tempreal=  -0 . 5+ (0 . 5*i/nprime) ; 


fiM 
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tempimag=  0 . 5- (0 . 5 *i/nprime) ; 

fprintf  (ofp,  "%i  %i  %e  %e\n",  a,  (npr:Lme-i) ,  tempreai,  tempimag) ; 

a=d+ 1 ; 

fcr  (.'c=0;  1<<  (nprime-i) ;  lc++)  {  /*  points  on  the  line  */ 

tempi=nprime-lc-l  ; 
temp]= (i+k) *sizes3/nprime+]; 
numreal=s3 [tempi] [tempj] .r; 
numiinag=s3  [tempi]  [temp]]  .i; 
fprintf (ofp, "%e  %e\n", numreal, numimag) ; 

}  /*  end  for  k=. . .  */ 

}  / *  end  for  ]= . . .  */ 

}  /*  end  for  i=. . .  */ 

/*  close  the  output  file  */ 
fclose (ofp) ; 

}  /*  end  if_else  */ 

}  /*  end  if  reduce=0...  */ 

if  (reduce==l)  {  /♦  reduce  the  amount  of  output  data  */ 

if  { ( (argc==8) && (argvl6] [!]==' a' ) ) I  I ( (argc==7) && (argv[5] [1]=='  a'  ) ) )  { 

/*  make  an  ascii  output  file  full  spectral  plane  */ 

/*  open  the  output  file  and  place  header  information  in  it  */ 
ofp  =  fopen (outf ile, "w" ) ; 
fprintf (ofp, "%i  %i\n", nprime,p) ; 
redindex=l; 

/*  write  out  data  values  to  file  for  plotting  */ 
for  (i=0;  i<nprime;  i++) { 
for  (j=»0;  j<p;  j++)  [ 

n;Knreal®s3 [i] [j*redindex] .r; 
numiinag«s3  [i]  [j*redindex]  .i; 

bigmag=sqrt (numreal *numreal  +  numimag*numimag) ; 
for  (k=l;  k<redindex;  k++)  [  /*  look  for  largest  value  */ 

tempreal=s3 [i] [ j*redindex  +k] .r; 
terapimag=s3 [i] [j*redindex  +k] .i; 

tempmag=sqrt  (tempreal*tenpreal  +  terrpimag*tenpimag) ; 
if  (tempmag>bigniag)  [  /*  found  a  bigger  value,  swap  */ 

b  iginag=t  empmag ; 
numreal=tempreal; 
numimag= t emp imag ; 

}  /*  end  if  tempmag>....  */ 

}  /*  end  for  k=l ...  */ 

fprintf (ofp, "%e  %e\n", numreal, numimag) ; 

}  /*  for  j=0 ...  */ 

}  /*  for  i=0 ...  */ 

/*  close  the  output  file  */ 
fclose (ofp) ; 

}  /*  if  ascii ...  */ 

else  [  /*  make  a  plot_sxaf  output  file  with  reduction  half  spectral  plane  */ 

/*  open  the  output  file  and  place  header  information  in  it  */ 
ofp  =  fopen (outfile, "w") ; 
redindex=sizes3/nprime; 
a*0;  /*  alpha  =  0  */ 

fprintf (ofp, "%i  %i  %e  %e  %i  %e  %e\n", 2, nprime, 0 . 0, 1 . 0, nprime, -0 . 5,  0 . 5) ; 
for  (i=0;  i<nprime;  i++)  [  /*  nprime  alpha  levels  */ 

tempreal-  -0 . 5+ (0 . 5*i/nprime) ; 
temp imag-  0 . 5- (0 . 5*i/nprime) ; 

fprintf (ofp, "%i  %i  %e  %e\n", a, (nprime-i) , tempreal, tempimag) ; 
a-a+1; 


as 
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for  (lc=0;  k<  (nprime-i)  ;  )<.++){  ■  *  points  at  the  alpha  level 

tempi=nprime-)c-l; 
temp]=(i+)c)  ’redindex; 
numreal=s3 [tempi J [temp]] .r; 
numimag=s3 [tempi] [temp]] .i; 

bigmag=sqrt (numreal*numreal  +  numimag*numimag) ; 
for  (]=1;  ]<redindex;  ]++)  {  /*  line  in  the  group  *■/ 

tempreal=s3 [tempi] [temp]  +]] .r; 
tempiraag=s3 [tempi] [temp]  +]] .i; 

tempmag=sqrt (tempreal*tempreai  +  tempimag'tempimag) ; 
if  (tempmag>bigmag)  [  /*  found  a  bigger  value,  swap  */ 

b  iginag=t  empmag  ; 
numreal=tempreal ; 
numimag=t  emp imag ; 

]  /*  end  if  tempmag>....  */ 

}  /*  end  for  ]=. . .  */ 

fprintf (ofp, "%e  %e\n", numreal, numimag)  ; 

}  /*  end  for  k=. . .  •/ 

}  /*  end  for  i=. . .  */ 

/*  close  the  output  file  */ 
f close (ofp) ; 

}  /*  end  if_else  */ 

}  /*  end  if  reduce=l...  */ 

} 


fib 


APPENDIX  E.  SUBFAM  PROGRAM  USE 


Correct  command  is; 

subfam  n  inputfilel  offsetl  freqzerol  conj 1  inputfile2  offsetC 
freqzero2  conj 2  outputfile  lOflag  Itypes 

where : 

n  is  the  number  of  samples  to  use  from  each  input file, 
it  must  be  a  power  of  2 

inputfilel  is  the  first  file  containing  signal  samples 

offsetl  is  the  initial  sample  position  offset  location 
in  inputfilel 

freqzerol  is  the  flag  indicating  which  frequency 
modification  to  perform  on  inputfilel,  options  are: 
zn  -  zero  neg  freqs,  shift  pos  freqs  down 
zp  -  zero  pos  freqs,  shift  neg  freqs  up 
zz  -  don't  do  anything  to  frequency  components 

conjl  is  the  flag  indicating  if  conjugation  of  the 
data  from  inputfilel  is  to  be  performed  prior  to 
multiplication  with  data  from  inputfile2,  options: 
y  -  yes,  conjugate  data 
n  -  do  not  conjugate  data 

inputfile2  is  the  second  file  containing  signal 
samples 

offset2  is  the  initial  sample  position  offset  location 
in  input file2 

freqzero2  is  the  flag  indicating  which  frequency 
modification  to  perform  on  inputfile2,  options  are: 
zn  -  zero  neg  freqs,  shift  pos  freqs  down 
zp  -  zero  pos  freqs,  shift  neg  freqs  up 
zz  -  don't  do  anything  to  frequency  components 

conj 2  is  the  flag  indicating  if  conjugation  of  the 
data  from  inputfile2  is  to  be  performed  prior  to 
multiplication  with  data  from  inputfilel,  options: 
y  -  yes,  conjugate  data 
n  -  do  not  conjugate  data 

outputfile  is  where  the  spectrum  values  will  be  placed 
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lOflag  indicates  the  datatypes  of  inputfilel, 

inputfile2  and  the  outputfile.  three  characters  are 
given  to  indicate  the  datatype  of  each  file, 
valid  datatypes  are: 
a  -  file  is  of  type  ascii 
b  -  file  is  of  type  binary 

ex.  "abb"  means  inputfilel  is  of  type  ascii, 
inputfile2  is  of  type  binary  and  the  outputfile 
is  to  be  of  type  binary 

Itypes  indicates  the  types  of  data  inside  inputfilel 
and  inputfile2.  the  outputfile  is  always  complex 
floating  point,  real  and  imaginary  components  on  the 
same  line.  two  letters  must  be  given  to  indicate 
the  datatype  for  each  input  file,  all  types  are 
converted  internally  to  complex  floating  point, 
valid  types  are: 

r  -  single  real  floating  point  samples  per  line 
i  -  single  integer  sample  per  line 
c  -  two  real  floating  point  numbers  per  line, 
representing  the  real  and  imaginary  parts 
of  a  single  sample 

ex.  "ic"  means  inputfilel  has  integer  samples 
and  inputfile2  has  complex  floating  point  samples 

example : 

subfam  8192  datal.dat  64  zn  n  data2.dat  2156  zp  y 
outputA.dat  bbb  rr 

input  file  format:  output  file  format: 

sample  #1  magnitude  #1  phase  #1 

sample  #2  magnitude  #2  phase  #2 


sample  #n  magnitude  #n  phase  #n 
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include  <stdlib.h> 

^include  <stdio.h> 

# include  <math.h> 

^include  " /homeS/carter/thesis/SSPI/pam/fft . c” 

^ include  " /home3 / cart er /thesis /SSP I /pam/radix.c" 

/*  values  used  in  up/down  conversion  */ 

^define  INTFREQ  49375000.0  /*  intermediate  signal  freq  in  Hz  ’/ 

^define  SAMPTIME  0.0000051282  /*  orig  signal  sampling  period  in  sec  ’/ 

main (argc, argv) 
int  argc; 
char  *argv[]; 

{ 

COMPLEX  *sl,*s2,*s3; 

int  1,  n, direction, norm,  of f 1, of f2,  quarter_n, scaninteger,  *readinteger,  *z; 
float  tempreal,  tempimag, outputmag, outputphase, numreal,  intfreq,  t,  convfac; 
float  *readreal, *readimag, *y; 

char  *infilel,  *infile2, *outfile, *con3l, *conj2, *freqzerol, *freqzero2,  *ioflag,  *itypes; 
FILE  *ifpl, *ifp2, *ofp; 

/* 

subfam  is  a  program  which  uses  signal  samples  from  two  input  files  to  calculate 
spectral  values  in  a  single  diamond  of  the  frequency-alpha  plane 

the  command  is  : 

subfam  n  inputfilel  offsetl  freqzerol  con^l  inputfile2  offset2  freqzero2 
conj2  outputfile  lOflag  Itypes 

n  is  the  number  of  samples  to  use  from  each  inputfile,  must  be  pwr  of  2 

inputfilel  is  the  name  of  the  first  input  file  containing  signal  samples 
offsetl  is  the  initial  sample  position  offset  location  in  inputfilel 
freqzerol  is  the  flag  indicating  which  frequency  modification  to  perform 
on  inputfilel,  valid  options  are: 

zn  -  zero  negative  frequency  components,  shift  positive 
coitponents  down 

zp  -  zero  positive  frequency  components,  shift  negative 
components  up 

zz  -  don't  do  anything  to  frequency  components 
con^l  is  the  flag  indicating  if  conjugation  of  the  data  from  inputfilel  is 
to  be  performed  prior  to  multiplication  with  data  from  inputfile2. 
valid  options  are: 
y  -  yes  conjugate  data 
n  -  do  not  conjugate  data 

inputfile2  is  the  name  of  the  second  input  file  containing  signal  saimples 
offset2  is  the  initial  sample  position  offset  location  in  inputfile2 
freqzero2  is  the  flag  indicating  which  frequency  modification  to  perform 
on  inputfile2,  valid  options  are: 

zn  -  zero  negative  frequency  components,  shift  positive 
components  down 

zp  -  zero  positive  frequency  components,  shift  negative 
components  up 

zz  -  don't  do  anything  to  frequency  components 
conj2  is  the  flag  indicating  if  conjugation  of  the  data  from  inputfile2  is 
to  be  performed  prior  to  multiplication  with  data  from  inputfilel. 
valid  options  are: 
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y  -  yes  conjugate  data 
n  -  do  not  con;]ugate  data 

outputfile  IS  the  filename  where  the  spectriim  values  will  be  placed 

lOflag  indicates  the  datatypes  of  inputfilel,  inputfile2  and  outputfile. 
three  letters  must  be  given  to  indicate  the  datatypes  for  all  files, 
valid  datatypes  are: 
a  -  file  is  of  type  ascii 
b  -  file  is  of  type  binary 
ex.  "abb"  means  inputfilel  is  of  type  ascii, 
input file2  is  of  type  binary,  and 
the  outputfile  is  to  be  of  type  binary 

Itypes  indicates  the  number  types  of  data  inside  inputfilel  and 

inputfile2.  the  outputfile  is  always  complex  floating  point, 
real  and  imaginary  components  on  the  same  line,  two  letters 
must  be  given  to  indicate  the  number  type  for  each  input  file, 
all  number  types  are  converted  to  complex  floating  point, 
valid  number  types  are: 

r  -  single  real  floating  point  sample  per  line 
i  -  single  integer  saitple  per  line 

c  -  two  floating  point  numbers  per  sample  on  each  line 
ex.  "ic"  means  inputfilel  has  integer  samples  and 

inputfile2  has  conplex  floating  point  samples 


sample  useage: 

subfam  8192  datal.dat  64  zn  n  data2.dat  2156  zp  y  outputA.dat  bbb  rr 
maintenance  history 

. . .  created  01  October  1992,  LCDR  Nancy  J.  Carter,  USN  at  NPGS  Monterey  CA 

...  10.15.92  -  added  freqzero  options  -  njc 

...  10.16.92  -  added  error  check  on  conjl  and  conj2  -  njc 

...  10.25.92  -  changed  down/up  shifting  to  same  code  as  in  old  version 

-  */ 

/*  check  for  the  correct  number  of  input  arguments  */ 

if  (argc  !=  13) { 

printf ("subfam. .. fatal  error. . .incorrect  number  of  calling  argumentsXn" ) ; 
printf ("\n") ; 

printf ("... .correct  format  is:  \n"); 

printf ("  subfam  n  inputfilel  offsetl  freqzerol  conjl  inputfile2  offset2  \n"); 


printf ("  freqzero2  conj2  outputfile  lOflag  ItypesXn"); 

printf ("\n") ; 
printf ("\n") ; 

printf  (" . n  is  the  number  of  samples  to  use  from  each  inputfile,  must  be  pwr  of  2\n") 

printf (" . inputfilel  is  a  file  containing  the  signal  samplesXn"); 

printf (" . offsetl  is  the  initial  sample  position  offset  in  inputfilelXn") ; 

printf (" . freqzerol  indicates  type  of  freq  zeroing  and  shifting  to  do  on  inputfilelXn 

printf ("  zn....zero  negative  freqs  and  shift  downXn"); 

printf  ("  zp....zero  positive  freqs  and  shift  upXn"); 

printf ("  z2....do  not  zero  anything,  do  not  shiftXn"); 

printf (" . conjl  indicates  if  conjugation  of  data  from  inputfilel  isXn"); 

printf  ("  desired  before  multiplicationXn")  ; 

printf  ("  y. . • .yes  conjugate  dataXn") ; 
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printfC  n....co  do  r.ot  conjugate  data\n"); 

printfC  usual  setup  is  ccn;]l  =  n,  001132  =  y\n''); 

print f(" . inputfile2  is  a  file  containing  the  signal  samples \n"); 

printfC . cff3et2  is  the  initial  sample  position  offset  in  inputf ile2\n'' )  ; 

printfC . freqzero2  indicates  type  of  freq  zeroing  and  shifting  to  do  on  inputfile 

printfC  zn....zero  negative  freqs  and  shift  downin''); 

printfC  zp....zero  positive  freqs  and  shift  up\n"); 

printfC  zz....do  not  zero  anything,  do  not  shiftin''); 

printfC . con]2  indicates  if  conjugation  of  data  from  inputfile2  is\n"); 

printfC  desired  before  multiplicationin" )  ; 

printfC  y....yes  conjugate  datain"); 

printfC  n....no  do  not  conjugate  datain"); 

printfC  usual  setup  is  conjl  =  n,  conj2  =  yin"); 

printfC . outputfile  is  where  the  spectrum  values  will  be  placedin"); 

printfC . lOflag  indicates  datatypes  of  inputfilel, inputf ile2  and  outputfilein" ) ; 

printf  ("  aaa . ascii, ascii,  asciiin") ; 

print  f  ("  aab . ascii,  ascii,  binary  in  " )  ; 

printf  ("  aba . ascii, binary,  asciiin" )  ; 

printf  ("  abb . ascii, binary, binaryin" ) ; 

printf  ("  baa . ascii, ascii, asciiin" )  ; 

printf  ("  bab . ascii,  ascii,  binaryin" )  ; 

printf  ( "  bba . binary,  binary,  asciiin" )  ; 

printf  ( "  bbb . binary,  binary,  binaryin" )  ; 

printfC . Itypes  indicates  number  types  of  input  samplesin"); 

printfC  r . single  real  floatin’’); 

printfC  i . single  integerin")  ; 

printfC  c . complex,  two  floatsin’’); 

printf C in") ; 

print  f  C . . . . sample  useage : in" ) ; 

printfC  subfam  8192  datal.dat  64  zn  n  data2.dat  2156  zp  y  outputA.dat  bbb  rrin’’); 
exit  0  ; 

} 

n=atoi (argv[l] ) ; 
inf ilel=argv [2 ]  ; 
offl=atoi (argv[3I )  ; 
f reqzerol=argv [ 4 ]  ; 
conjl=argv[5] ; 
inf ile2=argv [ 6] ; 
of f2=atoi (argv [7 ] ) ; 
f reqzero2=argv [ 8 ] ; 
con j2=argv [ 9] ; 
outfile=argv[10] ; 
ioflag=argv[ll] ; 
itypes=argv [ 12 ] ; 

/*  verify  that  n  is  a  power  of  2  */ 

i=n%2 ; 

If  (i'=0){ 

printf ( "subfam. . . fatal  errorin" ) ; 

printfC . calling  argument  n  is  not  a  power  of  2 in’’); 

exit  0  ; 

} 

/*  verify  that  conjl  &  conj2  are  y  or  n  •/ 

if  ( ( (conjl [0] !='y' )&& (conjl [0] !='n' ) )  I  I  ( (conj2 [0] !=' y' )&s (conj2 [0] !='n'  ) ) )  { 
printf ("subfam. . . fatal  errorin") ; 

printfC . conjl  and  conj2  must  be  set  to  y  or  nin’’); 


'la 
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exit  ( )  ; 

} 

/*  verify  that  itypes  are  r,  i  or  c  */ 

if  ( ( (itypes [0] ! =' r' ) && (itypes [0] ! =' i' ) && (itypes [0] ! ='  c'  ) )  I  I 
( (itypes [II ! =' r' ) && (itypes [1 ] !='!')&& (itypes [1 ] ! ='  c' ) ) )  { 
printf ("subfam. . .fatal  errorXn"); 


printfC . Itypes  must  be  set  to  r,  i  or  c\n''); 

exit  ( )  ; 

} 

/*  */ 

/*  allocate  space  for  the  1-D  data  arrays  */ 

/*  */ 


Sl= (COMPLEX*) calloc (n, sizeof (COMPLEX) )  ; 
if  (sl==NULL) { 

printf ("subfam.  .. insufficient  space  to  allocate  sl\n"); 
exit  0  ; 

} 

s2= (COMPLEX*) calloc (n, sizeof (COMPLEX) ) ; 

If  (s2==NULL) [ 

printf  ("subfcim.  ..  insufficient  space  to  allocate  s2\n"); 
exit  () ; 

} 

/*  allocate  space  to  read  in  binary  data  */ 

y= (float *)malloc (l*sizeof (float) ) ; 
z= (int*)malloc (l*sizeof (int) ) ; 

/*  */ 

/*  INPUT  */ 

/*  ♦/ 

/*  open  input  files,  get  the  signal  samples  */ 

/*  ♦/ 

switch  (ioflag[0])  { 

case  'a':  /*  inputfilel  is  an  ASCII  file  */ 

if  ((  ifpl  =  fopen (infilel, "r") )  ==  NULL  )  [ 

printf ("subfam. . .unable  to  open  ascii  inputfilelXn") ; 
exit  0  ; 

} 

/*  move  to  the  desired  starting  offset  position  in  the  file  */ 
for  (i=0; i<offl; i++) { 
switch  (itypes [0])  { 

case  ' r' :  /*  single  real  float  per  line  */ 

if  (fscanf  (ifpl,  "%e",  sterr^jreal)  !=1)  ( 

printf ("subfam. . .EOF  encountered  while  attempting  to  reach  offsetlXn"); 
fclose (ifpl) ; 
exit ( ) ; 

I 

break; 

case  ' i' :  /*  single  integer  per  line  */ 

if  (fscanf  (ifpl,  "%i",  scaninteger)  '.=!)  ( 

printf ("subfam. . .EOF  encountered  while  attenpting  to  reach  offsetlXn"); 
fclose (ifpl) ; 
exit  0  ; 

} 

break; 

default:  /*  coit^Jlex,  two  floats  per  line  */ 

if  (fscanf (ifpl, "%e  %e", &tempreal,&tempimag) !-2) [ 

printf ("subfam.  . .EOF  encountered  while  attempting  to  reach  offsetlXn"); 


13 
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fclose (ifpl) ; 
exit  0  ; 

} 

}  /’  end  switch  itypes[0]  '/ 

}  i*  end  if  i=0;i<offl  */ 

/*  read  in  n  data  samples  */ 
for  (1=0; i<n; 1++) { 
switch  (itypes[0])  ( 

case  ' r' :  /*  single  real  float  per  line  */ 

if  (fscanf (ifpl, "%e", &sl [i] .r) !=1) ( 

printf ("subfam.  . .EOF  reached  while  reading  inputfilelXn" ) ; 
fclose (ifpl) ; 
exit  0  ; 

) 

break; 

case  ' i' :  /*  single  integer  per  line  */ 

if  (fscanf (ifpl, "%i", scaninteger) !=1) { 

printf ("subfam. . .EOF  reached  while  reading  inputfilelXn”); 
fclose (ifpl) ; 
exit  0  ; 

} 

si [i] . r=scaninteger; 
break; 

default:  /*  complex,  two  floats  per  line  */ 

if  (fscanf (ifpl, "%e  %e", &sl [i] .r, isl [i] . i) 1=2) { 
printf ("subfam. . .EOF  reached  while  reading  inputfilelXn"); 
fclose (ifpl) ; 
exit  0 ; 

} 

}  /*  end  switch  itypes[0]  */ 

}  /*  end  for  i=0  i<n  •/ 

fclose (ifpl) ; 
break; 

case  'b' :  /*  inputfilel  is  a  BINARY  file  */ 

if  ((  ifpl  =  fopen (infilel, "rb") )  ==  NULL  )  { 

printf ( "subfam. . .unable  to  open  binary  inputfilelXn"); 
exit  ( )  ; 

} 

/♦  move  to  the  desired  starting  offset  position  in  the  file  */ 
for  (i=0; i<of f 1; i++) { 
switch  (itypes(0])  { 

case  ' r' :  /*  single  real  float  per  line  */ 

if  (fread(y, sizeof (*y) , 1, ifpl) !=1) { 

printf ("subfam. . .EOF  encountered  while  attempting  to  reach  offsetlXn"); 
fclose (ifpl) ; 
exit  0  ; 

) 

break; 

case  ' i' :  /*  single  integer  per  line  •/ 

if  (fread(z, sizeof (*z) , 1, ifpl) 1=1) ( 

printf ("subfam. . .EOF  encountered  while  attempting  to  reach  offsetlXn"); 
fclose (ifpl) ; 
exit  ( ) ; 
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break; 

default:  /*  complex,  two  floats  per  line  */ 

if  ( fread (y,  sizeof ( *y ) , 1, ifpl) ! =1) { 

print f ( "subfam. . .EOF  encountered  while  attempting  to  reach  offsetlXn"); 
fclose (ifpl) ; 
exit ( ) ; 

} 

if  (freacKy, sizeof (*y) , 1,  ifpl) (=1) { 

printf ("subfcun. . .EOF  encountered  while  attempting  to  reach  of fsetl\n" ) ; 
fclose (ifpl) ; 
exit  ( ) ; 

} 

}  /*  end  switch  itypes[0]  */ 

}  /*  end  if  i=0;i<offl  */ 

/*  read  in  n  data  samples  */ 
for  (i=0; i<n; i++) { 
switch  (itypes[0])  { 

case  ' r' :  /*  single  real  float  per  line  */ 

if  ( fread (y,  sizeof (*y) , 1, ifpl) !=1) ( 

printf  ("subfam.  .  .EOF  reached  while  reading  i.nputfilelVn"); 
fclose (ifpl) ; 
exit  0  ; 

si [i] .r=*y; 
break; 

case  ' i' :  /*  single  integer  per  line  */ 

if  ( fread (z, sizeof (*z) , 1, ifpl) ! =1)  1 

printf ("subfam. . .EOF  reached  while  reading  inputfilelXn") ; 
fclose (ifpl) ; 
exit  0  ; 

} 

si [i] . r=*z; 
break; 

default:  /*  complex,  two  floats  per  line  *■/ 

if  (fread(y,  sizeof (*y) , 1, ifpl) 1=1) { 

printf ("subfam.  . .EOF  reached  while  reading  inputf ilel\n" ) ; 
fclose (ifpl) ; 
exit  0  ; 

} 

si [i] .r=*y; 

if  (fread(y, sizeof (*y) , 1, ifpl) 1=1) { 

printf ("subfam. . .EOF  reached  while  reading  inputf ilel\n" ) ; 
fclose (ifpl) ; 
exit  0  ; 

} 

si [i] .r-*y; 

}  /*  end  switch  itypes[0]  */ 

)  /*  end  for  i=0  i<n  */ 

fclose  (ifpl) ; 
break; 

default:  /*  invalid  format  specified  for  inputfilel  */ 

printf ("subfaun. .. invalid  format  specified  for  inputf ilel\n" ) ; 
exit  () ; 

f  /♦  end  switch  ioflag[0]  */ 
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switch  (ioflag[l])  { 

case  ' a' :  /*  inputfile2  is  an  ASCII  file  */ 

if  ((  ifp2  =  fopen  (inf ile2, "r” ) )  ==  MULL  )  { 

printf ("subfam.  . .unable  to  open  ascii  inputf ile2\n" ) ; 
exit  ( )  ; 

} 

/*  move  to  the  desired  starting  offset  position  in  the  file 
for  (i=0; i<of f2 ; 1++) { 
switch  (itypes[l])  ( 

case  ' r' :  /•  single  real  float  per  line  */ 

if  (fscanf (ifp2, "he" , itempreai) !=1) { 

printf ("subfam. . .EOF  encountered  while  attempting  to  reach  of fset2\n" ) ; 
fclose (ifp2) ; 
exit  ( ) ; 

} 

break; 

case  ' i' :  /*  single  integer  per  line  */ 

if  (fscanf (ifp2, "%i", scaninteger) !=1) ( 

printf ("subfam. . .EOF  encountered  while  attempting  to  reach  oifset2\n") ; 
fclose (ifp2) ; 
exit  0  ; 

} 

break; 

default:  /*  complex,  two  floats  per  line  */ 

if  (fscanf  (ifp2,  "%e  %e",  iterttpreal,  stempimag)  !=2)  { 

printf ("subfam. . .EOF  encountered  while  attempting  to  reach  offset2\n"); 
fclose (ifp2) ; 
exit  0  ; 

} 

}  /♦  end  switch  itypes[l]  */ 

}  /*  end  if  i=0;i<off2  */ 

/*  read  in  n  data  samples  */ 
for  (i=0; i<n; 1++) { 
switch  (itypes[l])  ( 

case  ' r' :  /»  single  real  float  per  line  */ 

if  (fscanf (ifp2, "%e", &s2 [i] . r) ! =1) { 

printf ("subfam. . .EOF  reached  while  reading  inputf ile2\n" ) ; 
fclose (ifp2) ; 
exit  ( ) ; 

} 

break; 

case  '  i' :  /*  single  integer  per  line  */ 

if  (fscanf (ifp2, "hi" , scaninteger) !=1) { 

printf ("subfam.  . .EOF  reached  while  reading  inputf ile2\n" ) ; 
fclose (ifp2) ; 
exit  0  ; 

} 

s2 [i] . r=scaninteger; 
break; 

default:  /*  complex,  two  floats  per  line  */ 

if  (fscanf (ifp2, "%e  %e", &s2 [i] .r, &s2 ( i] . i) ! =2) ( 

printf ("subfam. . .EOF  reached  while  reading  inputfile2\n"); 
fclose (ifp2) ; 
exit  ( ) ; 


Dec  9  12:23  1992  testalg/subfam. c  Page  8 


}  /*  end  switch  itypes[lj  ” f 

}  /*  end  for  i=0  i<n  */ 

fclose (ifp2) ; 
break; 

case  'b';  /*  inputfile2  is  a  BINARY  file  */ 

if  ((  ifp2  =  f open (inf ile2, "rb") )  ==  NULL  )  { 

printf ("subfam. . .unable  to  open  binary  inputfile2\n" ) ; 
exit  ( )  ; 

} 

/♦  move  to  the  desired  starting  offset  position  in  the  file  */ 
for  (i=0; i<off2; i++) { 
switch  (itypes[l])  { 

case  ' r' :  /*  single  real  float  per  line  */ 

if  (fread(readreal, sizeof (*readreal) , 1, ifp2) !=1) { 

printf ("subfam. . .EOF  encountered  while  attempting  to  reach  offset2\n") 
fclose (ifp2 ) ; 
exit  0  ; 

} 

break; 

case  ' i' :  /*  single  integer  per  line  */ 

if  (fread(readinteger, sizeof (*readinteger) ,  1,  ifp2) !=1) { 

printf ("subfam. . .EOF  encountered  while  attenpting  to  reach  offset2\n") 
fclose (ifp2) ; 
exit ( ) ; 

break; 

default:  /*  complex,  two  floats  per  line  */ 

if  (fread(readreal, sizeof (*readreal) , 1, ifp2) !=1) { 

printf ("subfam. . .EOF  encountered  while  attempting  to  reach  offset2\n") 
fclose (ifp2) ; 
exit  0  ; 

} 

if  (fread(readimag, sizeof (*readimag) , 1, ifp2) !=1) { 

printf ("subfam. . .EOF  encountered  while  attempting  to  reach  offset2\n") 
fclose  (ifp2) ; 
exit  0  ; 

} 

}  /*  end  switch  itypes[l]  */ 

}  /♦  end  if  i=0;i<off2  */ 

/*  read  in  n  data  samples  ♦/ 
for  (i=0; i<n; i++) { 
switch  (itypes[l])  { 

case  ' r' :  /*  single  real  float  per  line  */ 

if  (fread(readreal, sizeof (*readreal) ,  1,  ifp2) !=1) { 

printf ("subfcun. . .EOF  reached  while  reading  inputfile2\n"); 
fclose (ifp2) ; 
exit  0  ; 

} 

s2[i] .r»*readreal; 
break; 

case  ' i' :  /*  single  integer  per  line  */ 

if  (fread(readinteger, sizeof (*readinteger) , 1, ifp2) !=1) { 

printf ("subfam. . .EOF  reached  while  reading  inputfile2\n"); 
fclose (ifp2) ; 
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exit  ( ) ; 

} 

scaninteger=*readinteger; 
s2  [  1  ]  .  r=scanir.teger; 
break.; 

default:  /*  complex,  two  floats  per  line  */ 

if  ( fread (readreal, sizeof ( *readreal) ,  1,  ifp2 )!=!){ 

printf ("subfam. . .EOF  reached  while  reading  inputf ile2\n" ) ; 
f close (ifp2) ; 
exit  0  ; 

} 

s2 [i] . r=*readreal; 

if  ( fread (readimag, sizeof (*readimag) ,  1,  ifp2) ! =1) { 

printf ("subfam. . .EOF  reached  while  reading  input f i le2 \ n ") ; 
fclose (ifp2) ; 
exit  ( ) ; 

s2[i] . i=*readimag; 

}  /*  end  switch  itypes[l]  */ 

}  /*  end  for  i=0  i<n  */ 

fclose (ifp2) ; 
break; 

default;  /♦  invalid  format  specified  for  inputfile2  */ 

printf ("subfam. .. invalid  format  specified  for  inputf ile2\n" ) ; 
exit  0  ; 

}  /*  end  switch  ioflag[l]  */ 

/*  */ 

/*  ZEROIZE  REQUESTED  FREQUENCY  SIDE  AND  SHIFT  APPROPRIATELY  ♦/ 

/*  */ 

quarter_n=n/4; 

switch  (freqzerol [1] )  { 

case  'n':  /*  zero  negative  and  shift  down  */ 

direction=l;  /*  request  a  forward  transform  time-to-freq  */ 
norm=0; 

f ft (si,  n,  direction,  norm) ;  /*  results  returned  in  si  */ 

for  (i=0;  i<(n/2);  i++) (  /*  zero  lower  half  */ 
si [i] .r=0; 
si [i] . i=0; 

} 

direction=-l;  /*  request  an  inverse  transform  freq-to-time  */ 
norm=2; 

f ft (si,  n,  direction,  norm) ;  /*  results  returned  in  si  */ 

/*  downconvert  */ 

t=0.0;  /*  starting  time  */ 

for  (i=0;  i<n;  i++)  { 

convf ac=TWOPI * INTFREQ*t ; 

tempreal=cos (convfac) ; 

tempimag= (-1) *sin (convfac)  ; 

numreal=sl [i] .r*tempreal  -  si [i] .i*tempimag; 

si  [i]  .  i=sl  [ i]  .  r*tentpiraag+sl  [i]  .  i*tempreal; 

sl[i] .r=numreal; 

t-t+SAMPTIME; 

} 

break; 

case  'p' :  /*  zero  positive  and  shift  up  */ 

direction*!;  /*  request  a  forward  transform  time-to-freq  */ 
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norin=0; 

f ft (si,  n,  direction, norm) ;  /*  results  returned  in  si  */ 

for  (i=(n/2);  i<n;  i++) {  /*  zero  upper  half  '! 
si [i] .r=0; 
si [i] .1=0; 

direction=-l;  /*  request  an  inverse  transform  freq-to-time  */ 
norm=2; 

f ft (si,  n, direction, norm) ;  /*  results  returned  in  si  */ 

/*  upconvert  */ 

t=0.0;  /*  starting  time  */ 

for  (i=0;  i<n;  i++)  { 

convfac=TWOPI*INTFREQ*t; 
terapreal=cos (convfac) ; 
tempimag=sin (convfac) ; 

numreal=sl [i] .r*tempreal  -  si [i] . i*tempimag; 
si [i] .i=sl [i] . r*tempimag+sl [i] .i*tempreal; 
sl[i] .r=numreal; 
t=t+SAMPTIME; 

} 

break; 

case  '  z' ;  /*  do  nothing  */ 

break; 

default:  /*  invalid  format  specified  for  freqzerol  */ 

printf ("subfam. .. invalid  format  specified  for  freqzerolVn" ) ; 
exit  ( )  ; 

}  /*  end  switch  */ 
switch  (freqzero2 [1] )  { 

case  'n':  /*  zero  negative  and  shift  down  */ 

direction=l;  /*  request  a  forward  transform  time-to-freq  */ 
norm=0; 

f ft (s2,  n,  direction, norm) ;  /*  results  returned  in  s2  */ 

for  (i=0;  i<(n/2);  i++) {  /*  zero  lower  half  */ 

s2 [i] .r=0; 
s2 [i] . i=0; 

} 

direction=-l;  /*  request  an  inverse  transform  freq-to-time  */ 
norm=2; 

f ft (s2,  n,  direction, norm) ;  /*  results  returned  in  s2  */ 

/♦  downconvert  */ 

t=0.0;  /*  starting  time  */ 

for  (i=0;  i<n;  i++)  { 

convf ac=TWOP I * INTFREQ*t ; 

tempreal=cos (convfac) ; 

tempimag= (-1) *sin (convfac) ; 

numreal=s2 [i] .r*tempreal  -  s2 [i] . i*tempimag; 

s2 [i] . i=s2 [i] .r*tempimag+s2 [i] .i*tempreal; 

s2[i] .r=numreal; 

t=t+SAMPTIME; 

1 

r 

break; 

case  'p' :  /*  zero  positive  and  shift  up  */ 

direction=l;  /*  request  a  forward  transform  time-to-freq  */ 
norm»0; 

f ft (s2, n, direction, norm) ;  /*  results  returned  in  s2  ♦/ 

for  (i-(n/2);  i<n;  i++) {  /*  zero  upper  half  */ 
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s2 [1]  .r=0; 

s2 [i: .1=0; 

} 

direccicn=-l ;  request  an  inverse  transform  freq-to-time  ’  / 

norm=2 ; 

f  ft  (s2,  n,  direction,  norm)  ;  /■'  results  returned  in  s2  ^ 
i'  upconvert 

t=0.0;  /*  starting  time 

for  (1=0;  i<n;  1++)  { 

convf ac=TWOP I  * INTFREQ*t ; 
tempreal=cos (convfac) ; 
tempimag=sin (convfac)  ; 

numreal=s2 [1] .r*tempreal  -  s2  [  1 ] . i^tempimag; 
s2 [ i] . i=s2 [ 1] . r*tempimag+s2 [1] . I’tempreal; 
s2 [1] . r=numreal; 
t  =t+ SAMP  TIME; 

} 

break; 

case  'z' :  /*  do  nothing  */ 

break; 

default:  /*  invalid  format  specified  for  freqzero2  */ 

printf ( "subfam. ,. invalid  format  specified  for  freqzero2\n" ) ; 
exit  ( )  ; 

}  /*  end  switch  */ 

/*  */ 

/*  MULTIPLICATION  */ 

/*  */ 

/*  create  space  for  the  output  array  */ 

/*  */ 

s3= (COMPLEX*) calloc (n, sizeof (COMPLEX) )  ; 
if  (s3==NULL) { 

printf ("subfam. .. insufficient  space  to  allocate  s3\n"); 
exit  ( )  ; 

} 

if  ( (con j 1 [0 ] ==' y' ) && (con]2 [0 ] ==' y' ) )  {  /*  multiply  conjugate  si  times  conjugate  s2  */ 
for  (i=0; i<n; 1++) i 

s3 [i] .r=sl [i] .r*s2 [i] .r-sl [i] .i*s2 [i] .1; 
s3[i] .i=(-l)*(sl[i] .i*s2[i] .i+sl[i] .i*s2[i] .r); 

} 

if  ( (con^l [ 0 ] ==' y' ) && (conj2 [0] ==' n' ) )  {  /*  multiply  conjugate  si  times  s2  */ 
for  (i=0; i<n; i++) { 

s3 [ i ] . r=sl [i] . r*s2 [i] . r+sl [1] .i*s2[i] .1; 
s3 [ i] . i=sl [i ] . r*s2 [ i] . i-sl [ i] . i*s2 [1] . r; 

} 

} 

if  ( (conjl [0] ==' n' ) && (conj2 [0] ==' y' ) )  {  /*  multiply  si  times  conjugate  s2  */ 
for  (i=0; i<n; i++) { 

s3 [i] .r=sl [i] .r*s2 [i] .r+sl [i] .i*s2 [1] .1; 
s3 [i] .i=sl  [i]  .i*s2  [i]  .r-sl [i] .r*s2 [i] .1; 

) 

> 

if  ( (con jl [0 ] ==' n' ) && (con j2 [0] ==' n' ) )  {  /*  multiply  si  times  s2  */ 
for  (i=0; i<n; i++) { 

s3 [i] .r=sl [i] .r*s2 [i] .r-sl [i] .i*s2[i] .i; 
s3 [i] .i=sl [i] .r*s2[i] .i+sl[i]  .i*s2[i] .r; 
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} 

} 

/*  take  an  FFT  of  n-point  s3  */ 

/-  »/ 

/*  set  values  in  arguments  sent  to  fft  */ 

direction=l;  /»  recjuest  a  forward  transform  time-to-freq  */ 

norm=0; 

:ft (s3,  n,  direction,  norm) ;  /*  results  returned  in  s3  */ 

*  / 

/*  OUTPUT  */ 

/*  */ 

/*  take  magnitude  and  phase  of  s3  and  place  in  the  output  file  * / 
I*  */ 

switch  (ioflag[2])  { 

case  ' a' :  /*  outputfile  is  an  ascii  file  */ 

if  ((  ofp  =  fopen (outf lie, "w" ) )  ==  NULL  )  { 

printf ("subfam. . .unable  to  open  ascii  outputf ile\n" ) ; 
exit  0  ; 

} 

/*  write  out  data  to  the  ascii  file  */ 

for  (i=0; i<n; i++) { 

outputmag=sqrt ((s3[i] .r*s3[i] .r)  +  (s3[i] .i*s3ti] .i)) ; 
outputphase=atan (s3 [i] . i/s3 [i] .r) ; 

fprintf (ofp, "%-16 . 6e  %-16 . 6e\n", outputmag, outputphase) ; 

} 

/*  close  the  output  file  */ 

fclose (ofp) ; 

break; 

case  'b' :  /*  outputfile  is  a  binary  file  */ 

if  ((  ofp  =  fopen (outfile, "wb") )  ==  NULL  )  ( 

printf ("subfam. . .unable  to  open  binary  outputfileXn") ; 
exit  0  ; 

/*  write  out  data  to  the  binary  file  */ 

for  (i=0; i<n; i++) ( 

outputmag=sqrt ((s3[i] .r*s3[i] .r)  +  (s3[i] .i*s3(i]  .i)); 
outputphase=atan (s3[i] .i/s3[i] .r); 

*y=outputraag; 

if  (fwrite (y, sizeof (*y) , l,ofp) !=1) { 

printf ("subfam. . .error  while  writing  outputf ile\n" ) ; 
fclose (ofp) ; 
exit  0  ; 

} 

*y=outputphase; 

if  (fwrite (y, sizeof (*y) , 1, ofp) ! =1) ( 

printf ("subfam. .. error  while  writing  outputfileXn" ) ; 
fclose (ofp) ; 
exit  0  ; 

} 

}  / *  end  if  i=0  ...  */ 

fclose (ofp) ; 

break; 

default:  /*  invalid  format  specified  for  outputfile  */ 

printf ("subfam. .. invalid  format  specified  for  outputfileXn") 
exit  ( ) ; 

}  /*  end  switch  ioflag[2]  */ 


LIST  OF  ABBREVIATIONS 


BPSK  binary  phase  shift  keying 

FAM  FFT  accumulation  method 

FFT  fast  fourier  transform 

FSM  frequency  smooching  method 

PAM  pulse  amplitude  modulation 

QAM  quadrature  amplitude  modulation 

QPSK  quadrature  phase  shift  keying 

SCF  spectral  correlation  function 

SSCA  strip  spectral  correlation  algorithm 

SSPI  Statistical  Signal  Processing,  Inc. 

SUBFAM  sub- FFT  accumulation  method 
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